2019-10-14 14:54:19 热度:

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

人全基因组分析一直是整个测序行业最重要的内容之一,随着各种测序仪性能的快速提升,人全基因组测序价格越来越便宜,周期越来越低,可以预见,即将有越来越多的人全基因组被测序出来。这里我们将系统介绍一个完整的人全基因组分析案例,人全基因组分析可以大致分为四个过程。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-anonymouscd hg38mirros 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_qcfastqc -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.htmlf

3、参考序列建立索引;

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

4、bwamem比对;

 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.samtime 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.bamsamtools 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.logtime 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.Rtime 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.Rtime 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.statsplot-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.bamsamtools view -h Sample1.sorted.bam | /ifs1/Software/biosoft/lumpy-sv-master/scripts/extractSplitReads_BwaMem -i stdin | samtools view -Sb - | samtools sort -> Sample1.splitters.sorted.bamlumpyexpress -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可视化变异结果。IntegrativeGenomicsViewer交互式基因组浏览器,它是一种高性能的可视化工具,用来交互式地探索大型综合基因组数据。它支持各种数据类型,包括array-based的和下一代测序的数据和基因注释。

上一篇:没有了

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

推荐阅读

基于金纳米粒子、外切酶放大信号、非靶点诱导的“行走”适配体开发的检测玉
玉米赤霉烯酮(Zearalenone,ZEN)存在于各类粮食作物如:大米和小麦,具有神经毒性、致癌性和免疫毒性,作为一种支原体雌激素会损害人们的健康。开发灵敏、准确、廉价的ZEN检测方法至关重要。基于适配体的比色传感器因其直观观察方便、简单、实时分析等优点,在常规分析中得到了广泛的应用。 作者利用核酸外切酶III(ExoIII)切割双链核苷酸的3-平末端或凹陷端的特性放大信号,再结合金纳米粒子的催化反应(4-硝基苯酚作为比色剂),并基于非靶点诱导的行走适配体开发了一种新型比色传感器,用于定量检测ZEN的含...[详细]
2019-11-11
一个人全基因组完整数据分析脚本
人全基因组分析一直是整个测序行业最重要的内容之一,随着各种测序仪性能的快速提升,人全基因组测序价格越来越便宜,周期越来越低,可以预见,即将有越来越多的人全基因组被测序出来。这里我们将系统介绍一个完整的人全基因组分析案例,人全基因组分析可以大致分为四个过程。1、从DNA到fastq;2、从fastq到bam;3、从bam到vcf;4、从vcf到pdf; 一、从DNA到fastq 从DNA到fastq也就是测序的过程,对于人全基因组的测序,要分清楚几个问题?1、选择哪种测序平台?2、需要提供多少样品?3、需要...[详细]
2019-10-14
  • 一个人全基因组完整数据分析脚本
    人全基因组分析一直是整个测序行业最重要的内容之一,随着各种测序仪性能的快速提升,人全基因组测序价格越来越便宜,周期越来越低,可以预见,即将有越来越多的人全基因组被测序出来。这里我们将系统介绍一个完整的人全基因组分析案例,人全基因组分析可以大致分为四个过程。1、从DNA到fastq;2、从fastq到bam;3、从bam到vcf;4、从vcf到pdf; 一、从DNA到fastq 从DNA到fastq也就是测序的过程,对于人全基因组的测序,要分清楚几个问题?1、选择哪种测序平台?2、需要提供多少样品?3、需要...
  • 基于金纳米粒子、外切酶放大信号、非靶点诱导的“行走”适配体开发的检测玉米赤霉烯酮比色传感器
    玉米赤霉烯酮(Zearalenone,ZEN)存在于各类粮食作物如:大米和小麦,具有神经毒性、致癌性和免疫毒性,作为一种支原体雌激素会损害人们的健康。开发灵敏、准确、廉价的ZEN检测方法至关重要。基于适配体的比色传感器因其直观观察方便、简单、实时分析等优点,在常规分析中得到了广泛的应用。 作者利用核酸外切酶III(ExoIII)切割双链核苷酸的3-平末端或凹陷端的特性放大信号,再结合金纳米粒子的催化反应(4-硝基苯酚作为比色剂),并基于非靶点诱导的行走适配体开发了一种新型比色传感器,用于定量检测ZEN的含...