全基因组完整数据实战

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

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

全基因组完整数据实战相关推荐

SPARK全栈全流程大数据实战之技术选型篇
2019独角兽企业重金招聘Python工程师标准>>> ###一.技术选型,环境搭建安装及生产部署 ####1.大数据研发调研和需求分析如果你已经或正在尝试搭建一套大数据环境或生 ...
linux中主成分分析软件,基于全基因组snp数据进行主成分分析（PCA）
现将如何基于全基因组的SNP数据进行PCA分析流程记录下来: 1)全基因组snp数据格式为 .vcf 2)利用vcftools软件进行格式转换(Linux系统下:进入 /vcftools_0.1.13 ...
基于全基因组测序数据鉴定结构变异的四大类方法总结
不同类型的基因组变异示意图(图片来源:labspaces) 上次给大家总结介绍了基因组单核苷酸多态性(single nucleotide polymorphism,SNP)的鉴定方法,今天给大家介绍结 ...
一个全基因组重测序分析实战
Original 2017-06-08 曾健明生信技能树这里选取的是 GATK best practice 是目前认可度最高的全基因组重测序分析流程,尤其适用于人类研究. PS:其实本文应该属于 ...
GWAS | 全基因组关联分析 | PLINK | 实战 | 统计遗传学
参考:PLINK | File format reference vcftools plink的主要功能:数据处理,质量控制的基本统计,群体分层分析,单位点的基本关联分析,家系数据的传递不平衡检验,多 ...
7分钟分析人类全基因组，他们刷新全球纪录，此前最快也要24小时
金磊发自凹非寺量子位报道 | 公众号 QbitAI 7分钟,这是来自中国的一支团队"合力出成绩".一举打破的世界纪录: 全球首次将人类全基因组分析,推进分钟级时代. 这支团 ...
GWAS | 原理和流程 | 全基因组关联分析 | Linkage disequilibrium (LD)连锁不平衡 | 曼哈顿图 Manhattan_plot | QQ_plot |...
问题: linkage disequilibrium (LD)和 pairwise correlation的区别?似乎它们都能达到相同的目的. 先从直觉上理解一下GWAS的原理: 核心就是SNP与表型 ...
【3月30日直播】新冠病毒全基因组测序——Midnight试剂盒及整体解决方案
识别上方二维码或点击「阅读原文」免费报名参加新冠疫情肆虐全球,基于Nanopore测序技术和数据分析在全球感染性疾病防控中的优势充分显现出来.该平台使用灵活.操作简便.产出快速.分析实时等特征为 ...
【生信】全基因组关联分析（GWAS）原理
[生信]全基因组关联分析(GWAS)原理文章的文字/图片/代码部分/全部来源网络或学术论文,文章会持续修缮更新,仅供大家学习使用. 目录 [生信]全基因组关联分析(GWAS) 1.前提知识介绍 1. ...

全基因组完整数据实战

全基因组完整数据实战相关推荐

最新文章

热门文章