人全基因组分析

基于基因学苑公众号的推送内容，将生信相关内容收集整理，便于后续查阅，仅做个人使用。

人全基因组分析可以大致分为四个过程。

从 DNA 到 fastq；
从 fastq 到 bam；
从 bam 到 vcf；
从 vcf 到 pdf；

从 DNA 到 fastq

从 DNA 到 fastq 也就是测序的过程，对于人全基因组的测序，要分清楚几个问题？
1、选择哪种测序平台？ 2、需要提供多少样品？ 3、需要多少测序量？如果想进行个人的全基因组测序，可以抽取 10ml 左右静脉血液即可；人的基因组大小是 3G 数据量，按照当前测序片段长度，例如双末端 150bp，需要测序 30 倍数据，也就是 90G 数据。

# 相关分析软件安装
#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

# 下载人类基因组数据库
lftp ftp.broadinstitute.org/bundle -u gsapubftp-anonymous
cd hg38
mirros hg38

从 fastq 到 bam

将测序的到的 fastq 文件，比对到参考序列，得到 bam 格式文件，需要对 bam 进行很多步处理。经过对 bam 一系列的处理，最终得到了经过排序，去除 duplication 以及 bqsr 之后的 bam 文件。

# 原始数据质控
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
# 数据过滤
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
# 参考序列构建索引
bwa index -p Homo_sapiens_assembly38 -a bwtsw Homo_sapiens_assembly38.fasta
gatk CreateSequenceDictionary -R Homo_sapiens_assembly38.fasta -O Homo_sapiens_assembly38.dict
# bam 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
# 标记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
# 碱基较正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

使用 gatk 进行变异检测，输入文件为 Sample1.sorted.markdup.BQSR.bam。

# 利用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
# 利用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 \
# 利用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
# 统计结果
bcftools stats Sample1.HC.snps.indel.VQSR.vcf.gz >  view.stats
plot-vcfstats view.stats -p output

其余突变检测

除了利用 gatk 找 SNP 和小的 InDel 之外，还可以利用 lumpy 找 SV 突变，CNVnator 找 CNV 突变。

# 利用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
# 利用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
# 利用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 文档报告。

# 利用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
# 利用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/
# 与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 可视化变异结果。
输入文件：

参考基因组
bam 文件
snp vcf 文件
indel vcf 文件

一个人全基因组完整分析流程

从 DNA 到 fastq

从 fastq 到 bam

从 bam 到 vcf

其余突变检测

从 vcf 到 pdf

可视化

CATALOG

FEATURED TAGS