欲练神功，必先得到VCF—变异信息提取（call SNP pip

群体遗传大部分的分析大部分基于VCF文件，所以得到一个高质量的VCF文件很有必要。在此记录一下从重测序的fastq文件提取SNP，过滤VCF的一系列流程。
之前的做法还是一步一步慢慢来，但涉及到更多的个体数据就很麻烦，因此写了一个小小的sh脚本，也算是节约时间吧。
废话不多说，上sh~
call_snp.sh 每一个小步骤需要改的地方不多，内存和线程根据服务器配置和需求改就行

###使用须知

#1）安装samtools bwa fastp gatk bedtools  （conda大法好）

#2）参考基因组只需建立一次索引 所有索引及ref在ref文件夹中

#3）更改ref等信息

ref=PATH_TO_REF

q1=$name_1.fq.gz

q2=$name_2.fq.gz

n=sample_name

#为参考基因组建立索引

bwa index $ref

samtools faidx $ref

gatk CreateSequenceDictionary -R $ref -O $ref.dict

#fastp 过滤

fastp -l 45 -q 20 -w 15 -i $q1 -I $q2 -o $n_clean_1.fq.gz -O $n_clean_2.fq.gz

#比对

bwa mem -t 30 -R "@RG\tID:$x\tPL:illumina\tLB:$x\tSM:$x" $ref $n_clean_1.fq.gz $n_clean_2.fq.gz | samtools view -Sb - > $n.bam

#排序

samtools sort -@ 30 -m 2G -O bam -o $n.sort.bam $n.bam

#标记重复序列

gatk MarkDuplicates -I $n.sort.bam  -O $n.sort.markdup.bam -M $n.sort.markdup_metrics.txt

samtools index $n.sort.markdup.bam

#变异检测

gatk HaplotypeCaller --java-options -Xmx30G  -R $ref --emit-ref-confidence GVCF -I $n.sort.markdup.bam -O $n.g.vcf

#分染色体运行 every_chr_gvcf.py  (一次运行代码如上，分染色体运行用every_chr_gvcf.py 生成的文件来算，优点：快就完事了)

#计算比对率 覆盖度 测序深度

samtools flagstat $n.sort.bam > $n.map

cat $n.map | head -n 5 | tail -n 1 | awk -F " " '{print$5}' | awk -F "(" '{print$2}' > out.map

bedtools genomecov -ibam $n.sort.bam -d > $n.depth

awk '{x+=$3} END {print x/NR}' $n.depth > $n.depth.1

awk '{print$3}' $n.depth | grep -w "0" | wc -l > $n.c.1

awk '{print$3}' $n.depth | wc -l > $n.c.2

paste $n.c.1 $n.c.2 > $n.c.3

awk '{print$1"\t"$2"\t"($2-$1)/$2}' 03  | awk '{print $3*100}' | sed 's/\(.\+\)/\1%/g' > $n.coverage

echo $n > name

paste -d "\t" out.map $n.depth.1 $n.coverage > $n.stat

chmod 777 $n.stat   ###最终统计文件是绿色的，内容是比对率 深度 覆盖度

every_chr_gvcf.py (脚本小垃圾 凑活看 能用但难看)

#make by JJ  python v3.7

import sys

print("usage: name.py chr.txt out sample_name path_to_ref"+"\n"+"mkdir chr")

f1 = open(sys.argv[1],'r')

f2 = open(sys.argv[2],'w')

n = sys.argv[3]

r = sys.argv[4]

n1 = str(n)

r1 = str(r)

l = f1.read().splitlines()

for i in range(len(l)):

        f2.write("gatk HaplotypeCaller --java-options -Xmx30G  -R "+r1+" --emit-ref-confidence GVCF -I "+n1+".sort.markdup.bam"+" -L "+str(l[i])+" -O chr/"+n1+str(l[i])+".g.vcf\n")

f1.close()

f2.close()

用到的软件：

Samtools

bwa

fastp

gatk 强烈推荐v 4.0

bedtools

后续步骤，且听下回。