2019-01-23_WES_实战总结

昨天晚上跑的 7E5240_raw.gvcf 跑完了，今天把丹丹跑的7E5239_bqsr.bam_raw.gvcf 7E5241_bqsr.bam_raw.gvcf复制到我的文件夹下了，之后就开始做下一步的合并gvf。

seq 22|while read id;do echo chr${id};done >bed.txt
cat >>bed.txt
chrX
chrY
chrM

先把所有染色体名称写进命名为 bed.txt的文档。
接下来执行合并样本的命令 ：

ref="/home/qmcui/database/reference/index/hisat/hg38/hg38.fa"
snp="/home/qmcui/dbsnp_146.hg38.vcf.gz"
cat bed.txt|while read bed;do gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" GenomicsDBImport -L $bed -R $ref -V 7E5239_bqsr.bam_raw.gvcf -V 7E5240_raw.gvcf -V 7E5241_bqsr.bam_raw.gvcf --genomicsdb-workspace-path gvcfs_${bed}.db
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" GenotypeGVCFs -R $ref -V gendb://gvcfs_${bed}.db --dbsnp $snp -O final_${bed}.vcf;done

此时出现第一报错：
-R 参数有问题
后来发现是忘了复制过来gatk的子命令GenomicsDBImport，-L/-R/-V都是这个子命令的参数，没有子命令这些参数自然就没有意义啦。
修改之后又出现第二个报错。

image.png
只是因为复制丹丹文件的时候只把7E5239_bqsr.bam_raw.gvcf 7E5241_bqsr.bam_raw.gvcf
复制到我的文件夹下了，其实每一个ID生成了3个文件，我们都应该复制过来，再运行就没问题了。
上一步运行完成之后生成了final_chr1.vcf
final_chr2.vcf
两个vcf文档。
其实正常情况应该都可以chr1-22/X/Y都可以生成，只是这里我们上课只生成两个意思一下。

### 合并vcf
gatk MergeVcfs -I final_chr2.vcf -I final_chr1.vcf -O raw.conbine.vcf.gz

过滤

mkdir vcf_filter
### vcf过滤
raw_vcf="/home/jhuang/wes/raw.conbine.vcf.gz"
filter_argu_snp="QUAL < 30.0 || QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR >
3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0"
gatk SelectVariants -select-type SNP -V $raw_vcf -O raw.snp.vcf.gz

gatk VariantFiltration -V raw.snp.vcf.gz --filter-expression $filter_argu_snp --filter-name "Filter" -O filter.snp.vcf.gz

gatk SelectVariants --exclude-filtered true -V filter.snp.vcf.gz -O filtered.snp.vcf.gz

上步过滤代码运行后会发生很诡异的事情，连报错都很不正常

image.png

后来丹丹交给我，不提前赋值，把所有需要赋值的地方都代入，再运行果然正常了。

下面接着进行过滤得到indel文件。

raw_vcf="/home/jhuang/wes/raw.conbine.vcf.gz"
filter_argu_snp="QUAL < 30.0 || QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR >
3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0"
gatk SelectVariants -select-type INDEL -V $raw_vcf -O raw.indel.vcf.gz

gatk VariantFiltration -V raw.indel.vcf.gz --filter-expression $filter_argu_indel --filter-name "Filter" -O filter.indel.vcf.gz

gatk SelectVariants --exclude-filtered true -V filter.indel.vcf.gz -O filtered.indel.vcf.gz

这步跟上一步一样，赋值的地方要直接代入，否则报错

 ####合并，也可不合并，看最后用途
gatk MergeVcfs -I filtered.snp.vcf.gz -I filtered.indel.vcf.gz -O conbine.filtered.vcf.gz

注释annovar

sample="7E5241"
/home/qmcui/software/ANNOVAR/annovar/convert2annovar.pl -format vcf4 -allsample filtered.indel.vcf.gz -outfile filtered.indel.vcf.gz.anno

/home/qmcui/software/ANNOVAR/annovar/annotate_variation.pl -buildver hg38 --outfile ${sample}.anno ${sample}.annovar /home/qmcui/software/ANNOVAR/annovar/humandb