一个人全基因组完整数据分析脚本

编者按
人全基因组分析一直是整个测序行业最重要的内容之一，随着各种测序仪性能的快速提升，人全基因组测序价格越来越便宜，周期越来越低，可以预见，即将有越来越多的人全基因组被测序出来。这里我们将系统介绍一个完整的人全基因组分析案例，人全基因组分析可以大致分为四个过程。
1、从DNA到fastq；
2、从fastq到bam；
3、从bam到vcf；
4、从vcf到pdf；

一、从DNA到fastq

从DNA到fastq也就是测序的过程，对于人全基因组的测序，要分清楚几个问题？
1、选择哪种测序平台？
2、需要提供多少样品？
3、需要多少测序量？

目前illumina的hiseq平台和bgiseq都可以完成人的全基因组测序，三代测序依然成本较大；如果想进行个人的全基因组测序，可以抽取10ml左右静脉血液即可；人的基因组大小是3G数据量，按照当前测序片段长度，例如双末端150bp，需要测序30倍数据，也就是90G数据；这里我们有一对完整的人全基因组测序原始数据；

-rwxr-xr-x 1 wangtong wangtong  24G 9月  19 22:51 180586B_4607-DNA_S33_L003_R1_001.fastq.gz
-rwxr-xr-x 1 wangtong wangtong  25G 9月  19 23:29 180586B_4607-DNA_S33_L003_R2_001.fastq.gz

二、分析准备工作

得到数据之后就可以开始进行分析；不过需要先进行一些准备工作，例如下载软件和数据库；
1、下载软件

#bwa:
https://jaist.dl.sourceforge.net/project/bio-bwa/bwa-0.7.17.tar.bz2
#samtools:
https://github.com/samtools/samtools/releases/download/1.8/samtools-1.8.tar.bz2
#GATK4:
https://github.com/broadinstitute/gatk/releases/download/4.0.2.1/gatk-4.0.2.1.zip
#bcftools:
https://github.com/samtools/bcftools/releases/download/1.8/bcftools-1.8.tar.bz2
#SNPeff:
https://jaist.dl.sourceforge.net/project/snpeff/snpEff_latest_core.zip
#lumpuy:
https://github.com/arq5x/lumpy-sv
#cnvnator:
https://github.com/abyzovlab/CNVnator/releases

2、下载数据库

lftp ftp.broadinstitute.org/bundle -u gsapubftp-anonymous
cd hg38
mirros hg38 

下载完数据，放到hg38/目录下；

├── 1000G_omni2.5.hg38.vcf.gz
├── 1000G_omni2.5.hg38.vcf.gz.tbi
├── 1000G_phase1.snps.high_confidence.hg38.vcf.gz
├── 1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi
├── Axiom_Exome_Plus.genotypes.all_populations.poly.hg38.vcf.gz
├── Axiom_Exome_Plus.genotypes.all_populations.poly.hg38.vcf.gz.tbi
├── dbsnp_146.hg38.vcf.gz
├── dbsnp_146.hg38.vcf.gz.tbi
├── hapmap_3.3_grch38_pop_stratified_af.vcf.gz
├── hapmap_3.3_grch38_pop_stratified_af.vcf.gz.tbi
├── hapmap_3.3.hg38.vcf.gz
├── hapmap_3.3.hg38.vcf.gz.tbi
├── Homo_sapiens_assembly38.fasta
├── Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
├── Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi
└── wgs_calling_regions.hg38.interval_list

三、从fastq到bam

测序数据，软件，数据库有了之后就可以开始进行分析了。接下来的工作就是将测序的到的fastq文件，比对到参考序列，得到bam格式文件，需要对bam进行很多步处理。
1、原始数据质控；

mkdir rawdata_qc
fastqc -f fastq -o rawdata_qc 180586B_4607-DNA_S33_L003_R1_001.fastq.gz  180586B_4607-DNA_S33_L003_R2_001.fastq.gz

2、数据过滤；

fastp -i 180586B_4607-DNA_S33_L003_R1_001.fastq.gz  -I 180586B_4607-DNA_S33_L003_R2_001.fastq.gz -o wgs_clean.1.fq.gz  -O wgs_clean.2.fq.gz  -z 4 -q 20 -u 30 -n 6 -w 4 -h clean.html
f

3、参考序列建立索引；

bwa index -p Homo_sapiens_assembly38 -a bwtsw Homo_sapiens_assembly38.fasta
gatk CreateSequenceDictionary -R Homo_sapiens_assembly38.fasta -O Homo_sapiens_assembly38.dict

4、bwa mem比对；

 bwa mem -t 4 -M -Y -R "@RG\tID:Sample1\tPL:ILLUMINA\tLB:Lib1\tSM:Sample1" hg38/Homo_sapiens_assembly38 wgs_clean.1.fq.gz wgs_clean.2.fq.gz >Sample1.sam

time gatk SortSam -I Sample1.sam -O Sample1.sorted.bam -SO coordinate --CREATE_INDEX true

#也可以利用samtools进行排序
#samtools sort -@ 4 -o Sample1.sorted.bam Sample1.sam
#rm -rf Sample1.sam

5、标记duplication；

gatk MarkDuplicates -I Sample1.sorted.bam -M Sample1.markdup_metrics.txt -O Sample1.sorted.markdup.bam
samtools index Sample1.sorted.markdup.bam
#rm -rf Sample1.sorted.bam

6、碱基较正BQSR

#加time命令计时
time gatk BaseRecalibrator \
         -R hg38/Homo_sapiens_assembly38.fasta \
         -I Sample1.sorted.markdup.bam \
         --known-sites hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
         --known-sites hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
         --known-sites hg38/dbsnp_146.hg38.vcf.gz \
         -O Sample1.sorted.markdup.recal_data.table >bqsr.log

time gatk ApplyBQSR \
         --bqsr-recal-file Sample1.sorted.markdup.recal_data.table \
         -R hg38/Homo_sapiens_assembly38.fasta \
         -I Sample1.sorted.markdup.bam \
         -O Sample1.sorted.markdup.BQSR.bam 
time samtools index Sample1.sorted.markdup.BQSR.bam

#rm -rf Sample1.sorted.markdup.bam

四、从bam到vcf

经过对bam一系列的处理，最终得到了经过排序，去除duplication以及bqsr之后的bam文件，接下来就可以使用gatk进行变异检测，输入文件为Sample1.sorted.markdup.BQSR.bam；
1、利用gatk得到vcf文件

time gatk HaplotypeCaller \
   --emit-ref-confidence GVCF \
   -R hg38/Homo_sapiens_assembly38.fasta \
   -I Sample1.sorted.markdup.BQSR.bam \
   -O Sample1.HC.g.vcf.gz 
time gatk GenotypeGVCFs \
   -R hg38/Homo_sapiens_assembly38.fasta \
   -V Sample1.HC.g.vcf.gz \
   -O Sample1.HC.vcf.gz 

2、利用VQSR方法过滤SNP结果

time  gatk VariantRecalibrator \
            -R hg38/Homo_sapiens_assembly38.fasta \
            -V Sample1.HC.vcf.gz \
            --resource hapmap,known=false,training=true,truth=true,prior=15.0:hg38/hapmap_3.3.hg38.vcf.gz \
            --resource omni,known=false,training=true,truth=false,prior=12.0:hg38/1000G_omni2.5.hg38.vcf.gz \
            --resource 1000G,known=false,training=true,truth=false,prior=10.0:hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
            --resource dbsnp,known=true,training=false,truth=false,prior=2.0:hg38/dbsnp_146.hg38.vcf.gz \
            -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP \
            -mode SNP \
            -O Sample1.HC.snps.recal \
            --tranches-file Sample1.HC.snps.tranches \
            --rscript-file Sample1.HC.snps.plots.R

time gatk ApplyVQSR \
         -R hg38/Homo_sapiens_assembly38.fasta \
         -V Sample1.HC.vcf.gz \
         -O Sample1.HC.snps.VQSR.vcf.gz \
         --recal-file Sample1.HC.snps.recal \
         --tranches-file Sample1.HC.snps.tranches \
         -mode SNP \

3、利用VQSR处理InDel

time gatk VariantRecalibrator \
             -R hg38/Homo_sapiens_assembly38.fasta \
             -V Sample1.HC.snps.VQSR.vcf.gz \
             --max-gaussians 4 \
             --resource mills,known=false,training=true,truth=true,prior=12.0:hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
             --resource dbsnp,known=true,training=false,truth=false,prior=2.0:hg38/dbsnp_146.hg38.vcf.gz \
             -an QD -an DP -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
             -mode INDEL \
             -O Sample1.HC.snps.indel.recal \
             --tranches-file Sample1.HC.snps.indel.tranches \
             --rscript-file Sample1.HC.snps.indel.plots.R


time gatk ApplyVQSR \
                 -R hg38/Homo_sapiens_assembly38.fasta \
                 -V Sample1.HC.snps.VQSR.vcf.gz \
                 -O Sample1.HC.snps.indel.VQSR.vcf.gz \
                 --truth-sensitivity-filter-level 99.0 \
                 --tranches-file Sample1.HC.snps.indel.tranches \
                 --recal-file Sample1.HC.snps.indel.recal \
                 -mode INDEL

最终，得到的Sample1.HC.snps.indel.VQSR.vcf.gz既是需要的文件。

4、统计结果

#统计突变数目
bcftools stats Sample1.HC.snps.indel.VQSR.vcf.gz >  view.stats
plot-vcfstats view.stats -p output

五、其余突变检测

除了利用gatk找SNP和小的InDel之外，还可以利用lumpy找SV突变，CNVnator找CNV突变；
1、利用lumpy检测SV突变

samtools view -b -F 1294 Sample1.sorted.bam  | samtools sort - > Sample1.discordants.sorted.bam
samtools view -h Sample1.sorted.bam | /ifs1/Software/biosoft/lumpy-sv-master/scripts/extractSplitReads_BwaMem -i stdin | samtools view -Sb - | samtools sort -> Sample1.splitters.sorted.bam
lumpyexpress -B Sample1.sorted.bam -S Sample1.discordants.sorted.bam -D Sample1.splitters.sorted.bam -o Sample1.lumpu.sv.vcf

2、利用delly检测SV突变

#SV检测
delly call -g hg38/Homo_sapiens_assembly38.fasta -o Sample1.delly.sv.bcf -n Sample1.sorted.bam
#过滤结果
delly filter -f germline -p -q 20 Sample1.delly.sv.bcf -o Sample1.delly.sv.filter.bcf

3、利用CNVnator检测CNV突变

#1.提取mapping信息
/ifs1/Software/biosoft/CNVnator/cnvnator -root Sample1.root -tree Sample1.sorted.bam -unique 

#2.生成质量分布图HISTOGRAM
/ifs1/Software/biosoft/CNVnator/cnvnator -root Sample1.root -his 100  -d hg38/Homo_sapiens_assembly38.fasta 

#3.生成统计结果
/ifs1/Software/biosoft/CNVnator/cnvnator -root Sample1.root -stat 100 

#4.RD信息分割partipition
/ifs1/Software/biosoft/CNVnator/cnvnator -root Sample1.root -partition 100 

#5.变异检出
/ifs1/Software/biosoft/CNVnator/cnvnator -root Sample1.root -call 100 > Sample1.cnvnator.vcf

六、从vcf到pdf

该阶段主要是对的到的突变vcf文件Sample1.HC.snps.indel.VQSR.vcf.gz,与各种数据库进行比对注释，得到注释结果之后，利用LaTex或者html语言进行标记，最终生成PDF文档报告。

1、利用SNPeff进行注释

#列出所有数据库    
java -jar /ifs1/Software/biosoft/snpEff/snpEff.jar databases | less     
#筛选人基因组数据库    
java -jar /ifs1/Software/biosoft/snpEff/snpEff.jar  databases |grep "Homo"  
java -jar  /ifs1/Software/biosoft/snpEff/snpEff.jar -i vcf -o vcf GRCh38.86 Sample1.HC.snps.indel.VQSR.vcf.gz >Sample1.snpeff.vcf

2、利用Annovar进行注释

#利用annovar进行注释

#1 装换格式
/ifs1/Software/biosoft/annovar/convert2annovar.pl -format vcf4old Sample1.HC.snps.indel.VQSR.vcf.gz >Sample1.annovar.input
#2 进行注释
/ifs1/Software/biosoft/annovar/annotate_variation.pl -buildver hg38 --geneanno --outfile Sample1.anno Sample1.annovar.input /ifs1/Software/biosoft/annovar/humandb/ 

3、与Clinvar数据库注释

#clinvar数据库注释
#perl annotate_variation.pl -downdb -webfrom annovar clinvar_20180603 -buildver hg38 humandb/
/ifs1/Software/biosoft/annovar/convert2annovar.pl -format vcf4old Sample1.HC.snps.indel.VQSR.vcf.gz >Sample1.annovar.input
/ifs1/Software/biosoft/annovar/annotate_variation.pl --filter -buildver hg38 --outfile Sample1.clinvar.anno Sample1.annovar.input -dbtype clinvar_20180603 /ifs1/Software/biosoft/annovar/humandb/

其他一些常用数据库：
HGMD：https://www.qiagenbioinformatics.com/hgmd-resources/
SNPedia：https://www.snpedia.com/
PhramGKB: https://www.pharmgkb.org/

七、突变结果可视化

很多时候，为了更加精细化的查看突变结果，可以利用IGV可视化变异结果。Integrative Genomics Viewer 交互式基因组浏览器，它是一种高性能的可视化工具，用来交互式地探索大型综合基因组数据。它支持各种数据类型，包括array-based的和下一代测序的数据和基因注释。

1、下载安装
下载地址：http://software.broadinstitute.org/software/igv/home
需要Java：
2、输入文件：
参考基因组
bam文件
snp vcf文件
indel vcf文件

我也可以分析人全基因组数据吗？

当然可以，我们已经将以上所有数据，软件，数据库，案例脚步放到服务器里，购买服务器账号之后，即可开始进行练习，体验一次完整的全基因组数据分析。

$ ls -1
a1.index.sh
a2.fasqc.sh
a3.fastp.sh
a4.bwa.sh
a5.markdup.sh
a6.bqsr.sh
b1.gatk.sh
b2.vqsr_snp.sh
b3.vqsr_indel.sh
b4.hardfilter.sh
b5.delly.sh
b6.lumpy.sh
b7.cnvnator.sh
b8.stat.sh
c1.snpEff.sh
c2.annovar.sh
c3.SnpSift.sh
c4.clinvar.sh

想要更加系统的学习全基因组数据怎么办？

可以参加我们11月28日-12月2日，举行的生物信息培训班（南京站），里面会系统介绍人全基因组分析，期待你的到来。