GATK4流程學(xué)習(xí)之背景知識(shí)與前期準(zhǔn)備 - 簡(jiǎn)書(shū)
GATK4流程學(xué)習(xí)之DNA-Seq variant calling(Germline：SNP+INDEL) - 簡(jiǎn)書(shū)
GATK4流程學(xué)習(xí)之RNA-Seq variant calling(SNP+INDEL) - 簡(jiǎn)書(shū)
補(bǔ)：Mutect2+scRNAseq+cancer cell - 簡(jiǎn)書(shū)
從四個(gè)樣本（human sample）的全基因組測(cè)序(WGS)結(jié)果fastq.gz開(kāi)始选脊，通過(guò)
GATK Best Practices系列分析叙身，得到基于參考基因組的樣本SNP+INDEL變異信息侥蒙，以VCF文件為最終結(jié)果洞辣。筆記中主要使用之前搭建的conda環(huán)境
conda activate GATK
mkdir dna-seq
cd dna-seq
GATK Best Practices
一、流程概括
- 1啰脚、下載原始測(cè)序文件fastq.gz彤灶。如果下載到的是sra格式，還需要進(jìn)行一步轉(zhuǎn)換操作席覆。
- 2、對(duì)fastq.gz進(jìn)行質(zhì)控汹买。如有必要娜睛，還要進(jìn)行必要的過(guò)濾。
- 3卦睹、測(cè)序信息比對(duì)到參考基因組，得到原始的bam格式文件方库。
- 4结序、對(duì)比對(duì)結(jié)果的bam文件進(jìn)行“優(yōu)化”。主要包括去重復(fù)與校準(zhǔn)堿基質(zhì)量分?jǐn)?shù)兩步纵潦。
- 5徐鹤、根據(jù)優(yōu)化好的多樣本bam結(jié)果進(jìn)行variant calling，最終得到VCF結(jié)果邀层。
二返敬、下載數(shù)據(jù)
示例數(shù)據(jù)來(lái)自D. melanogaster WGS paired-end Illumina data with NCBI accessions SRR1663608,SRR1663609, SRR1663610, SRR1663611, corresponding to samples ZW155, ZW177, ZW184, and ZW185 respectively.
法1、直接下載fastq.gz格式(推薦)
- 數(shù)據(jù)源來(lái)自EBI寥院，https://www.ebi.ac.uk/ena/browser/view/SRR1663608
- 使用軟件aspera-cli
mkdir raw_fq
cd raw_fq
cat > fq.txt
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR166/008/SRR1663608/SRR1663608_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR166/008/SRR1663608/SRR1663608_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR166/009/SRR1663609/SRR1663609_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR166/009/SRR1663609/SRR1663609_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR166/000/SRR1663610/SRR1663610_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR166/000/SRR1663610/SRR1663610_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR166/001/SRR1663611/SRR1663611_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR166/001/SRR1663611/SRR1663611_2.fastq.gz
cat fq.txt |while read id
do
ascp -QT -l 300m -P33001 \
-i ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh \
era-fasp@$id .
done
download result
- 由于只是學(xué)習(xí)操作劲赠，而一般全基因組測(cè)序結(jié)果很大，后續(xù)流程比較耗時(shí)。所以這里對(duì)每個(gè)樣本測(cè)序結(jié)果抽取10%序列進(jìn)行演示凛澎。
- 軟件：seqtk
mkdir ../raw_fq_sub
ls *fastq.gz | cut -d "." -f 1 > sample.txt
cat sample.txt | while read id
do
seqtk sample ${id}.fastq.gz 0.1 | gzip > ../raw_fq_sub/${id}_sub1.fastq.gz
done
ls lh
法2霹肝、下載sra格式，再轉(zhuǎn)為fastq.gz
- 數(shù)據(jù)源來(lái)自NCBI塑煎，https://www.ncbi.nlm.nih.gov/sra/?term=SRR1663608
- 使用軟件 sra-tools ,可直接根據(jù)提供的SRR號(hào)下載
mkdir raw_fq
cd raw_fq
cat fq.txt
prefetch SRR1663608
prefetch SRR1663609
prefetch SRR1663610
prefetch SRR1663611
#格式轉(zhuǎn)換
fastq-dump --split-files *.sra --gzip
試了
prefetch方法沫换，結(jié)果下載失敗了，不知道是不是因?yàn)榫W(wǎng)絡(luò)的原因最铁。還是推薦第一種方法讯赏。
image.png
三、原始測(cè)序數(shù)據(jù)質(zhì)控
1冷尉、質(zhì)控
- 軟件：fastq
fastqc -t 10 ./*sub1.fastq.gz -o ./qc.out
-
如果質(zhì)控結(jié)果不好漱挎，可進(jìn)行序列過(guò)濾。
image.png
- 如上圖結(jié)果网严，所有的測(cè)序結(jié)果基本都合格识樱。除了read后面的堿基測(cè)序質(zhì)量不咋樣，因此也可不進(jìn)行過(guò)濾震束。
- 如果需要過(guò)濾的話怜庸，軟件有很多，這里使用的是trimmomatic
mkdir ../clean_fq
ls *sub1.fastq.gz |cut -d"_" -f 1 | cut -d"/" -f 2 | uniq| while read id
do
trimmomatic PE -phred33 -threads 10 \
${id}_1_sub1.fastq.gz ${id}_2_sub1.fastq.gz \
../clean_fq/out.${id}_1_sub1.fastq.gz ../clean_fq/out.trim.${id}_1_sub1.fastq.gz \
../clean_fq/out.${id}_2_sub1.fastq.gz ../clean_fq/out.trim.${id}_2_sub1.fastq.gz \
ILLUMINACLIP:/home/shensuo/biosoft/Trimmomatic/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 \
SLIDINGWINDOW:4:5 LEADING:5 TRAILING:5 MINLEN:25
done
#out***文件為過(guò)濾結(jié)果
#out.trim***文件為被去除點(diǎn)的不合格序列信息
- 過(guò)濾參數(shù)釋義
（1）ILLUMINACLIP是為了刪除fastq里的接頭信息垢村。在手動(dòng)安裝trimmomatic軟件時(shí)割疾，可找到提供的參考接頭信息。
（2）SLIDINGWINDOW：當(dāng)檢測(cè)到某read序列出現(xiàn)連續(xù)4個(gè)堿基時(shí)嘉栓，把包括這4個(gè)堿基在內(nèi)的后面所有堿基切除掉宏榕。
（3）LEADING與TRAILING分別表示切除read首尾質(zhì)量低于5的堿基。
（4）MINLEN表示若read長(zhǎng)度小于25侵佃，則跳過(guò)麻昼。
根據(jù)trim的結(jié)果，可再次進(jìn)行fastqc檢測(cè)馋辈。由于原始數(shù)據(jù)的QC結(jié)果基本合格抚芦，故后續(xù)分析還是基于
raw_fq_sub文件夾的結(jié)果進(jìn)行分析。
raw_fq_sub
四迈螟、測(cè)序序列比對(duì)到參考基因組
- 軟件：bwa 【better for DNA-seq】
mkdir bam
#參考數(shù)據(jù)庫(kù)索引文件
#由于人類基因組較大叉抡，故建索引較耗時(shí)，大概需要1-2h答毫。因此在準(zhǔn)備工作時(shí)已建好
BWA_INDEX=/home/shensuo/biosoft/GATK/gatk-4.1.9.0/bundle/hg38/bwa_index/gatk_hg38
ls ./raw_fq_sub/*gz | cut -d"/" -f 3 | cut -d"_" -f 1 | uniq | while read id
do
bwa mem -t 30 -M -R "@RG\tID:${id}\tSM:${id}\tLB:${id}\tPL:Illumina" $BWA_INDEX \
./raw_fq_sub/${id}_1_sub1.fastq.gz ./raw_fq_sub/${id}_2_sub1.fastq.gz | \
samtools view -Sb - > ./bam/${id}.bam ;\
samtools sort -@ 10 -m 4G -O bam -o ./bam/${id}.sorted.bam ./bam/${id}.bam
done
#先利用bwa比對(duì)得到sam文件褥民，再利用samtools轉(zhuǎn)為bam二進(jìn)制文件
- bwa參數(shù)解釋
（0）memMaximal Exact Match：Reports only one alignment for each read
（1）-M表示如果一條read的不同部分有不同的最合適比對(duì)位點(diǎn)，則分別記錄各個(gè)part的最佳位點(diǎn)洗搂。
（2）-t表示多線程比對(duì)
（3）-R表示為測(cè)序結(jié)果補(bǔ)充read group信息@RG消返，再后續(xù)操作步驟中起到重要作用载弄，因此著重說(shuō)明下 -
@RG
（1）一般理想NGS測(cè)序過(guò)程是：一張測(cè)序玻片有8條lane，每條lane為1個(gè)樣本的1個(gè)文庫(kù)進(jìn)行測(cè)序侦副，每條lane的測(cè)序過(guò)程可以理解為獨(dú)立的侦锯。
image.png
所以ID就是每條lane的編號(hào)(一定確保unique)，SM就是樣本名(sample)秦驯，LB為構(gòu)建的文庫(kù)名尺碰，PL為測(cè)序平臺(tái)，目前主流都是illumina
@RG
（2）復(fù)雜的測(cè)序?qū)嶒?yàn)設(shè)計(jì)包括①同一sample建不同文庫(kù)译隘，多次測(cè)序亲桥。
情況1
②一個(gè)文庫(kù)包括多個(gè)樣本(序列所加標(biāo)簽不同)在一條lane里同時(shí)測(cè)序。如下圖是文庫(kù)里包括2個(gè)sample固耘，測(cè)了兩次情況下的@RG記錄情況
情況2
在示例代碼中题篷，統(tǒng)一用各自的SRR號(hào)代替。具體分析根據(jù)實(shí)際情況設(shè)置厅目。
關(guān)于sam/bam文件的意義及格式番枚，可參考之前的筆記http://www.reibang.com/p/f0f1f293f0bd
五、比對(duì)結(jié)果bam優(yōu)化之去重復(fù)
1损敷、why
image.png
-
由于Optical duplicates: (Illumina)葫笼、 Library duplicates (aka PCR duplicates)兩種誤差因素導(dǎo)致在建庫(kù)及測(cè)序時(shí)，導(dǎo)致一些reads含量“假陽(yáng)性”地增高拗馒，會(huì)影響后續(xù)variant calling的結(jié)果（False Positive variant call）
image.png - 因此這一步的目的就是通過(guò)一些方法路星，找到這些“異常表達(dá)”的read，并進(jìn)行去重處理至一個(gè)正常的水平诱桂。
2洋丐、how
- 軟件：GATK4的
MarkDuplicates
cd bam
ls *sorted.bam | cut -d"." -f 1-2 | while read id
do
gatk MarkDuplicates \
-I ${id}.bam \
-O ${id}.dedup.bam \
-M ${id}.dedup_metrics.txt ; \
samtools index ${id}.dedup.bam
done
rm *[0-9].bam
rm *[0-9].sorted.bam
其實(shí)GATK4的
MarkDuplicates功能是來(lái)自Picard軟件，在4.0版本時(shí)挥等，將功能集成到gatk里了友绝，便于用戶使用。在后面的步驟中肝劲，基本上用的都是GATK自己的功能了迁客。image.png
note
- 如上圖所示，此時(shí)應(yīng)該執(zhí)行INDEL Realignment涡相。主要目的是檢查是否有不合適的比對(duì)結(jié)果。
- 因?yàn)椴煌谋葘?duì)結(jié)果剩蟀，往往會(huì)出現(xiàn)不同的變異信息（如下圖）催蝗，而B(niǎo)WA只會(huì)從最優(yōu)比對(duì)的角度考慮，可能會(huì)沒(méi)有發(fā)現(xiàn)最有可能存在的變異信息育特。
- 這一步就是使用新的比對(duì)算法(GATK)丙号，并結(jié)合已知變異數(shù)據(jù)庫(kù)先朦，對(duì)BWA比對(duì)結(jié)果所產(chǎn)生的變異位置進(jìn)行局部重比對(duì)。
- Re-alignment no longer recommended if the genotyping method used downstream involves local haplotype assembly犬缨，例如后面使用的HaplotypeCaller喳魏。 所以這一步就省略了。
image.png
六怀薛、比對(duì)結(jié)果bam優(yōu)化之校準(zhǔn)堿基質(zhì)量分?jǐn)?shù)
1刺彩、why
- 堿基質(zhì)量分?jǐn)?shù)就是指測(cè)序fastq結(jié)果的堿基質(zhì)量值。
- 校正的前提是測(cè)序錯(cuò)誤堿基的可能性遠(yuǎn)遠(yuǎn)大于基因組變異的概率枝恋，或者說(shuō)物種基因變異很写淳蟆；把所有與參考基因組不一致的堿基視為測(cè)序錯(cuò)誤導(dǎo)致（可排除已知的變異記錄）焚碌。
- 然后根據(jù)得到的堿基質(zhì)量分?jǐn)?shù)畦攘，進(jìn)行轉(zhuǎn)換得到堿基數(shù)目，比如說(shuō)100個(gè)質(zhì)量分?jǐn)?shù)為20(10的-2次方)的堿基=1個(gè)錯(cuò)誤堿基十电。
- 最后根據(jù)一定的算法/模型知押，進(jìn)行質(zhì)量分?jǐn)?shù)的校正（如下圖），具體可深入了解鹃骂。
- 所以這里強(qiáng)調(diào)是堿基質(zhì)量分?jǐn)?shù)校準(zhǔn)這一步適合于變異概率很少台盯，并且有已知參考變異數(shù)據(jù)庫(kù)的物種基因組。因此這一步基本上只適合人類的測(cè)序數(shù)據(jù)偎漫。
image.png
測(cè)序的系統(tǒng)誤差也是校準(zhǔn)的重點(diǎn)之一爷恳。每條lane的測(cè)序過(guò)程是相對(duì)獨(dú)立的，即需要分別對(duì)來(lái)自不同lane的測(cè)序結(jié)果單獨(dú)校正象踊，這也是之前在生成sam文件温亲，需要添加
@RG的原因。
2杯矩、how
- 軟件：GATK4
- 分為兩步栈虚，首先
BaseRecalibrator利用已有的snp數(shù)據(jù)庫(kù)，建立相關(guān)性模型史隆，產(chǎn)生重校準(zhǔn)表( recalibration table)魂务，輸入已知的變異位點(diǎn)數(shù)據(jù)庫(kù)，用于屏蔽那些不需要重校準(zhǔn)的部分泌射。 - 然后
ApplyBQSR根據(jù)上一步得到的模型對(duì)原始?jí)A基質(zhì)量分?jǐn)?shù)進(jìn)行調(diào)整粘姜，且只會(huì)調(diào)整非已知SNP區(qū)域。
ls *sorted.dedup.bam | cut -d"." -f 1-3 | while read id
do
gatk BaseRecalibrator \
-I ${id}.bam \
-R /home/shensuo/biosoft/GATK/gatk-4.1.9.0/bundle/hg38/Homo_sapiens_assembly38.fasta \
--known-sites /home/shensuo/biosoft/GATK/gatk-4.1.9.0/bundle/hg38/dbsnp_146.hg38.vcf.gz \
--known-sites /home/shensuo/biosoft/GATK/gatk-4.1.9.0/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
--known-sites /home/shensuo/biosoft/GATK/gatk-4.1.9.0/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
-O ${id}.recal_data.table ; \
gatk ApplyBQSR \
--bqsr-recal-file ${id}.recal_data.table \
-R /home/shensuo/biosoft/GATK/gatk-4.1.9.0/bundle/hg38/Homo_sapiens_assembly38.fasta \
-I ${id}.bam \
-O ${id}.recal.bam
done
rm *dedup.bam*
六熔酷、Varint calling
-
在完成上述的處理后孤紧，就可以對(duì)每個(gè)樣本的比對(duì)結(jié)果(*.bam)提取變異信息了；
image.png -
如上圖是先利用pairHMM算法得到的每個(gè)樣本的gvcf變異記錄拒秘，然后再匯總所有樣本的信息号显，最后轉(zhuǎn)為VCF結(jié)果臭猜。
image.png
step1 對(duì)每個(gè)樣本產(chǎn)生 gvcf 文件(十分耗時(shí)，建議后臺(tái)運(yùn)行)
mkdir ../vcf
cd ../vcf
nohup ls ../bam/*.recal.bam | cut -d"/" -f 3 | cut -d"." -f 1 | while read id
do
gatk HaplotypeCaller \
--native-pair-hmm-threads 10 \
-R /home/shensuo/biosoft/GATK/gatk-4.1.9.0/bundle/hg38/Homo_sapiens_assembly38.fasta \
--dbsnp /home/shensuo/biosoft/GATK/gatk-4.1.9.0/bundle/hg38/dbsnp_146.hg38.vcf.gz \
-I ../bam/${id}.sorted.dedup.recal.bam \
--minimum-mapping-quality 30 \
-ERC GVCF \
-O ${id}.g.vcf
done > log.out 2>&1 &
#花費(fèi)了4.5h
--dbsnp參數(shù)表示如果有與參考變異數(shù)據(jù)庫(kù)一致的變異記錄押蚤，則將數(shù)據(jù)庫(kù)該記錄的ID編號(hào)整合到gvcf里（vcf格式的第2行）蔑歌。
image.png
step2 合并gvcf，combine g.vcf files
gatk CombineGVCFs \
-R /home/shensuo/biosoft/GATK/gatk-4.1.9.0/bundle/hg38/Homo_sapiens_assembly38.fasta \
--variant SRR1663608.g.vcf \
--variant SRR1663609.g.vcf \
--variant SRR1663610.g.vcf \
--variant SRR1663611.g.vcf \
-O SRR4.g.vcf
step3 轉(zhuǎn)為vcf揽碘，joint variant calling with GenotypeGVCFs
gatk GenotypeGVCFs \
-R /home/shensuo/biosoft/GATK/gatk-4.1.9.0/bundle/hg38/Homo_sapiens_assembly38.fasta \
-V SRR4.g.vcf \
-stand-call-conf 5 \
-O SRR4.vcf
七次屠、之后的處理(未代碼演示)
1、過(guò)濾
得到原始VCF文件后钾菊，首先要過(guò)濾到低質(zhì)量的Variant帅矗。有如下兩種方法
（1）硬過(guò)濾
- 根據(jù)VCF記錄的各個(gè)角度描述的變異評(píng)價(jià)信息指標(biāo)進(jìn)行過(guò)濾
- gatk的
VariantFiltration功能
image.png
（2）軟過(guò)濾
- Machine learning model, recommended instead of hard (threshold-based) filtering when a set of true, reliable variants is available。
- 當(dāng)?shù)玫降脑糣CF文件里包含一些與參考變異數(shù)據(jù)庫(kù)相同的calling時(shí)煞烫』氪耍可通過(guò)對(duì)這些variant的指標(biāo)情況建立模型，去評(píng)價(jià)raw vcf里未知的變異記錄滞详。
- gatk的
VariantRecalibrator,ApplyVQSR功能
image.png
2凛俱、探索
- 得到高質(zhì)量的vcf后，就可以基于VCF格式的理解料饥，對(duì)該結(jié)果進(jìn)行不同角度的探索蒲犬，比如2號(hào)染色體的特定位點(diǎn)的變異信息；質(zhì)量值高于100的變異信息......
-
vcftools軟件則提供了分析VCF結(jié)果的功能岸啡，可進(jìn)一步深入學(xué)習(xí)下原叮。http://vcftools.sourceforge.net/
image.png
