一、同时处理批量文件的话，先准备样本名称：
ls *R1* > 1
ls *R2* > 2
paste 1 2 > config
二、质控，去接头
cat config | while read id
do
arr=(${id})
fq1=${arr[0]}
fq2=${arr[1]}
echo $fq1 $fq2
# fastqc $fq1 $fq2
trim_galore --phred33 --fastqc --stringency 3 -q 20 --length 28 --paired $fq1 $fq2
done
——————————————————————————
mkdir Rawdata RawdataQC TrimmedFQ TrimmedQC TrimmedReport
将文件各自归类备用
三、FQ to VCF
用bwa将去过接头的数据比对到参考基因组上得到sam文件 | 得到bam文件 | 得到排序的bam文件
${path to bwa}bwa mem -t 16 ${path to reference}genome.fa $fq1 $fq2 | ${path to samtools}samtools view -@ 16 -bS | ${path to samtools}samtools sort -@ 16 > ${fq1%%_*}.sorted.bam
-t/@：指定线程数
对排序的bam文件建立索引
${path to samtools}samtools index ${fq1%%_*}.sorted.bam > ${fq1%%_*}.sorted.bam.bai
标记重复
picard MarkDuplicates INPUT=${fq1%%_*}.sorted.bam OUTPUT=${fq1%%_*}.sorted.dedup.bam METRICS_FILE=${fq1%%_*}.metrics.txt REMOVE_DUPLICATES=true
SNP calling
${path to samtools}samtools mpileup -Q 20 -q 20 -l ${your bed file} -D -ugf ${path to reference}genome.fa ${fq1%%_*}.sorted.dedup.bam | ${path to bcftools}bcftools call -m > ${fq1%%_*}.dedup.D.vcf
-D：输出单个样本深度
四、根据需要对vcf文件按照家系进行合并：
bgzip -c $var > $var.gz
tabix -f $var.gz
{path to your bcftools}bcftools merge -O v --force-samples $var1.vcf.gz $var2.vcf.gz > pedi.merge.vcf
网友评论