嘗試一下jimmy的WES流程旅东，雖然現(xiàn)在云生信勢頭正猛男娄，做生信的門檻對普通人來說降低了行贪，對專業(yè)人士提高了許多，但是誰能阻止我們學(xué)習(xí)的步伐呢模闲？相信那句「沒有白走的路」建瘫，繼續(xù)學(xué)習(xí)，學(xué)習(xí)是我的愛好和樂趣围橡。
這是生信技能樹上帖子的地址暖混，上面有所有代碼，贊翁授！一個(gè)如此熱愛分享的大神拣播！http://www.biotrainee.com/thread-122-1-1.html
一晾咪、軟件和數(shù)據(jù)準(zhǔn)備
1.軟件安裝
我的原則是能用conda解決的不去wget，所以就conda解決了大部分軟件安裝問題贮配。我的系統(tǒng)黑蘋果10.11.5（EI Capitan）谍倦，遠(yuǎn)景論壇上下載的大神配置好的clover配置文件，基本完美泪勒，我們的目的不是裝系統(tǒng)昼蛀，而是用嘛，所以也不深入研究了圆存。
thinkpad-E431
處理器：Intel Core i3-3120M 2*2.49 GHz
內(nèi)存：12 GB
conda install bwa
#bwa是基于BWT的生物序列比對工具:bwa是為二代測序的短序列比對參考序列（reference叼旋，fasta格式）而開發(fā)的比對軟件。
conda install samtools
#samtools是目前廣泛使用額關(guān)于處理bam/sam文件的工具沦辙。從測序儀得到的fastq文件經(jīng)過mapping后得到二進(jìn)制的bam文件夫植，而sam則是它的十進(jìn)制版本。
conda install gatk
#gatk主要用于從sequencing 數(shù)據(jù)中進(jìn)行variant calling油讯，包括SNP详民、INDEL。比如現(xiàn)在風(fēng)行的exome sequencing找variant陌兑，一般通過BWA+GATK的pipeline進(jìn)行數(shù)據(jù)分析沈跨。
wget https://software.broadinstitute.org/gatk/download/auth?package=GATK
#上面只是安裝了個(gè)gatk-register，仍需要下載gatk用gatk-register安裝
tar -jxvf GenomeAnalysisTK-3.8-0.tar.bz2
cd GenomeAnalysisTK-3.8-0-ge9d806836/
mv GenomeAnalysisTK.jar ../../biosofts/
gatk-register ~/biosofts/GenomeAnalysisTK.jar
#拷貝到目標(biāo)文件夾(安裝目錄)
conda install picard
#picard是一套操作高通量測序（HTS）數(shù)據(jù)的命令行工具(java)兔综，格式如SAM/BAM/CRAM和VCF饿凛。這些文件格式是在Hts規(guī)范存儲(chǔ)庫中定義的。尤其是SAM規(guī)范和VCF規(guī)范软驰。
2.數(shù)據(jù)準(zhǔn)備
1）模擬產(chǎn)生數(shù)據(jù)
因?yàn)閷?shí)際的數(shù)據(jù)太大了笤喳，對我這種弱弱的計(jì)算機(jī)簡直是一種折磨，我下過全部的sra序列碌宴，轉(zhuǎn)換成fastq后來比對，比對就是花了幾小時(shí)蒙畴，什么動(dòng)靜也沒有贰镣。所以，對于學(xué)習(xí)來說膳凝，產(chǎn)生一點(diǎn)模擬數(shù)據(jù)是不錯(cuò)的碑隆，jimmy人性化的為我們提供了pirs這個(gè)軟件來模擬產(chǎn)生測序數(shù)據(jù)。
我表示裝這個(gè)軟件折騰了良久蹬音，不過終于可以使用了上煤。首先這個(gè)軟件依賴于zlib和boost編譯，我才開始嘗試下載源碼編譯著淆，好像是成功了劫狠。但是拴疤，我不放心，于是用brew安裝了一下測試了一下独泞。
brew install zlib
brew install boost
然后是pirs這個(gè)軟件的安裝：
大概我只是按步驟裝了下就可以用了呐矾，中間有點(diǎn)報(bào)錯(cuò)，是sh文件里有個(gè)路徑錯(cuò)誤懦砂，我還以為是軟件的嚴(yán)重錯(cuò)誤蜒犯。
wget ftp://ftp.genomics.org.cn/pub/pIRS/pIRS_111.tgz #下載軟件
tar zxvf pIRS_111.tgz #解壓
cd pIRS_111/
make #提示說沒有源碼不需要編譯于是嘗試運(yùn)行了一下，竟然可以運(yùn)行荞膘。
mkdir -p data reference
pirs simulate -i reference/chr1.fa \
-s ~/biosofts/pirs/pIRS_111/Profiles/Base-Calling_Profiles/humNew.PE100.matrix.gz \
-d ~/biosofts/pirs/pIRS_111/Profiles/GC-depth_Profiles/humNew.gcdep_200.dat \
-b ~/biosofts/pirs/pIRS_111/Profiles/InDel_Profiles/phixv2.InDel.matrix \
-l 100 -x 50 -Q 33 -o A
pirs軟件參數(shù)學(xué)習(xí)：
simulate：模擬illumina reads數(shù)據(jù)
-i：參考基因組序列
-s <double> 雜合 SNP rate (default=0.001)
-d <double> GC內(nèi)容覆蓋深度罚随？
-b <string> input InDel-error profile,input InDel-error profile for simulating InDel-error of reads, (default=(exe_path)/Profiles/InDel_Profiles/phixv2.InDel.matrix)
-l:read讀長
-x:測序覆蓋深度，默認(rèn)5
-Q:33/64
-o <string> output prefix (default=Illumina)
#程序運(yùn)行過程：
Start to preview error profile...
In Base-calling profile:
Ref_Base_num: 4
Statistical_Cycle_num: 200
Seq_Base_num: 4
Quality_num: 41
100bp reads total average substitution error rate: 0.00372389
Start to construct simulation matrix...
Loading file: /Users/zd200572/biosofts/pirs/pIRS_111/Profiles/Base-Calling_Profiles/humNew.PE100.matrix.gz
Have finished constructing Base-calling simulation matrix1
Start to construct simulation matrix...
Loading file: /Users/zd200572/biosofts/pirs/pIRS_111/Profiles/Base-Calling_Profiles/humNew.PE100.matrix.gz
Have finished constructing Base-calling simulation matrix2
Loading the GC bias profile: /Users/zd200572/biosofts/pirs/pIRS_111/Profiles/GC-depth_Profiles/humNew.gcdep_200.dat
Have finished constructing GC bias simulation matrix
Start to construct InDel-error matrix...
Loading file: /Users/zd200572/biosofts/pirs/pIRS_111/Profiles/InDel_Profiles/phixv2.InDel.matrix
For simulating GC bias:
GC% abundance
0 0.00422544
1 0.00624158
2 0.00912628
...
In InDel-error profile:
profile_Cycle_num: 200
read1 count: 147091454
read2 count: 145785569
length of max InDel: 3
insertion-base rate of 100-bp read1: 4.77512e-06
deletion-base rate of 100-bp read1: 1.11823e-05
insertion-base rate of 100-bp read2: 8.36749e-06
deletion-base rate of 100-bp read2: 1.55069e-05
Have finished constructing InDel-error simulation matrix
Have finished reading scaffold chr1
Begin to output reads...
Output 10000 pair reads
Output 20000 pair reads
...
Finish output reads
#由于50x速度實(shí)在生成的太慢羽资，2x的試試淘菩，哈哈475s搞定。
In ouput data:
sum of Insertion in read1 file: 1211
sum of Deletion in read1 file: 2814
sum of Insertion in read2 file: 2103
sum of Deletion in read2 file: 3825
All done! Run time: 475s.
#對比一下50x的
In ouput data:
sum of Insertion in read1 file: 30374
sum of Deletion in read1 file: 70211
sum of Insertion in read2 file: 52633
sum of Deletion in read2 file: 95478
All done! Run time: 11248s.
2）真實(shí)的wes數(shù)據(jù)（不是必需）
沒看到在流程中用到削罩，下下來以后備用瞄勾，有100G之巨。測試了一下弥激，對于我的電腦配置进陡，單單bwa比對運(yùn)行要花100個(gè)小時(shí)。微服。趾疚。
# ~100GB disk space required to download the NA12878
wget [url]ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR262/ERR262997/ERR262997_1.fastq.gz[/url]
wget [url]ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR262/ERR262997/ERR262997_2.fastq.gz[/url]
3）參考基因組
cd reference
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.chromFa.tar.gz
tar zxvf hg38.chromFa.tar.gz
cd chroms
for i in $(seq 1 22) X Y M; do cat chr${i}.fa >> hg19.fasta; done
cp chr1.fa ../
rm -fr chr*.fasta
cd ../
4）基因組index和dictionary的獲得
bwa index -p hg38.chr1 chr1.fa
samtools faidx hg38.fa
picard CreateSequenceDictionary.jar R=hg38.fasta O=hg38.dict
5）其他參考數(shù)據(jù)
# dbSNP138
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/database/snp150.txt.gz
# GATK resource bundle
GATK_FTP=ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/2.8/hg38
# download library file for GATK's VariantRecalibrator
for f in hapmap_3.3.hg38.vcf.gz 1000G_omni2.5.hg38.vcf.gz 1000G_phase1.indels.hg38.vcf.gz Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
do
curl --remote-name ${GATK_FTP}/${f}
done
二、軟件運(yùn)行
1.bwa比對
由于對于一個(gè)軟件來說以蕴，功能太多糙麦，不可能一一覆蓋，于是丛肮，我就先學(xué)習(xí)這次用到的參數(shù)赡磅，順便對軟件做一個(gè)大概了解，由于是初學(xué)宝与，目標(biāo)先定低一點(diǎn)焚廊，先了解大概，以后做項(xiàng)目再一一學(xué)習(xí)习劫。
bwa的mem參數(shù)：MEM是bwa中最新的算法咆瘟，基本取代了前兩種，目的是將測序reads或組裝的contig比對到Reference上诽里。（Bio-Xin博客：http://www.cnblogs.com/leezx/p/6226717.html）
-t參數(shù)袒餐，線程數(shù)；-R參數(shù)：設(shè)置reads標(biāo)頭，“\t”分割灸眼；M——將較短的split hits標(biāo)記為secondary卧檐，與picard兼容；后邊跟參考基因組幢炸、reads文件和>以及要生成的sam文件泄隔。（如果GATK call SNP 必須用-r 參數(shù)）
R1=/mnt/A_100_500_1.fq.gz
R2=/mnt/A_100_500_2.fq.gz
#align with BWA, with 4 threads (-t 4). Reading Group is required (-R ...).
bwa mem -t 4 -R '@RG\tID:group_id\tPL:illumina\tSM:sample_id' $DATA/hg19.chr1/hg19.chr1 ${R1} ${R2} > A.chr1.sam
Real time: 457.713 sec; CPU: 1489.602 sec
#2x的數(shù)據(jù)
Real time: 12440.568 sec; CPU: 43493.514 sec #50x
2.picard
1）SortSam工具
由BWA生成的sam文件時(shí)按字典式排序法進(jìn)行的排序（lexicographically）進(jìn)行排序的（chr10，chr11…chr19宛徊，chr1佛嬉，chr20…chr22，chr2闸天，chr3…chrM暖呕，chrX，chrY）苞氮，但是GATK在進(jìn)行callsnp的時(shí)候是按照染色體組型（karyotypic）進(jìn)行的（chrM湾揽，chr1，chr2…chr22笼吟，chrX库物，chrY），因此要對原始sam文件進(jìn)行reorder贷帮∑萁遥可以使用picard-tools中的ReorderSam完成，同時(shí)格式轉(zhuǎn)換成了bam撵枢。（https://www.plob.org/article/7009.html）
picard SortSam CREATE_INDEX=True SORT_ORDER=coordinate VALIDATION_STRINGENCY=LENIENT INPUT=A.chr1.sam OUTPUT=A.chr1.sort.bam
#2x運(yùn)行情況
picard.sam.SortSam done. Elapsed time: 4.54 minutes.
Runtime.totalMemory()=872415232
#50x
2）MarkDuplicates工具
過量擴(kuò)增的reads并不是基因組自身固有序列民晒，不能作為變異檢測的證據(jù)，因此锄禽，要盡量去除這些由PCR擴(kuò)增所形成的duplicates潜必，這一步可以使用picard-tools來完成。去重復(fù)的過程是給這些序列設(shè)置一個(gè)flag以標(biāo)志它們沃但，方便GATK的識(shí)別磁滚。還可以設(shè)置 REMOVE_DUPLICATES=true 來丟棄duplicated序列。對于是否選擇標(biāo)記或者刪除宵晚，對結(jié)果應(yīng)該沒有什么影響恨旱，GATK官方流程里面給出的例子是僅做標(biāo)記不刪除。這里定義的重復(fù)序列是這樣的：如果兩條reads具有相同的長度而且比對到了基因組的同一位置坝疼，那么就認(rèn)為這樣的reads是由PCR擴(kuò)增而來，就會(huì)被GATK標(biāo)記谆沃。
picard MarkDuplicates CREATE_INDEX=true REMOVE_DUPLICATES=True ASSUME_SORTED=True VALIDATION_STRINGENCY=LENIENT METRICS_FILE=/dev/null INPUT=A.chr1.sort.bam OUTPUT=A.chr1.sort.dedup.bam
#2x運(yùn)行情況
picard.sam.markduplicates.MarkDuplicates done. Elapsed time: 5.15 minutes.
Runtime.totalMemory()=822083584
3）FixMateInformation
這一步驟在網(wǎng)上搜到了這個(gè)钝凶。。。（看來我得好好理出一個(gè)現(xiàn)在使用的流程來耕陷，生物信息的流程發(fā)展還真是快呀５嗝）在一些老版本的pipeline中(如seqanswers上的那個(gè))，這一步做完后還要對pair-end數(shù)據(jù)的mate information進(jìn)行fix哟沫。不過新版本已經(jīng)不需要這一步了饺蔑，Indel realign可以自己fix mate。（http://www.bbioo.com/lifesciences/40-115982-1.html）
picard FixMateInformation SO=coordinate VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true INPUT=A.chr1.sort.dedup.bam OUTPUT=A.chr1.sort.dedup.mate.bam
#2x運(yùn)行情況
picard.sam.FixMateInformation done. Elapsed time: 6.36 minutes.
Runtime.totalMemory()=877658112
3.gatk
終于到了大頭gatk出場了嗜诀，我表示我對它好陌生猾警，幾乎是第一次運(yùn)行它。通過上面可以看出他是一個(gè)經(jīng)常更新的程序隆敢，而且它雖然是一個(gè)java程序发皿，但是它有好幾個(gè)子工具，就像上面的picard拂蝎。先看看它的流程穴墅，官方的：
1)RealignerTargetCreator
對INDEL周圍進(jìn)行realignment，這個(gè)操作有兩個(gè)步驟温自，第一步生成需要進(jìn)行realignment的位置的信息, 輸出一個(gè)包含著possible indels的文件玄货。
GenomeAnalysisTK -R reference/chr1.fa \
-T RealignerTargetCreator \
-I A.chr1.sort.dedup.mate.bam -o A.intervals
2)IndelRealigner
第二步才是對這些位置進(jìn)行realignment, 最后生成一個(gè)realin好的BAM文件sample.realn.bam
GenomeAnalysisTK -R reference/chr1.fa \
-T IndelRealigner \
-targetIntervals A.intervals \
-I A.chr1.sort.dedup.mate.bam -o A.chr1.sort.dedup.mate.relaign.bam
3)BaseRecalibrator
（對base的quality score進(jìn)行校正）堿基質(zhì)量分?jǐn)?shù)重校準(zhǔn)（Base quality score recalibration，BQSR)悼泌，就是利用機(jī)器學(xué)習(xí)的方式調(diào)整原始堿基的質(zhì)量分?jǐn)?shù)松捉。它分為兩個(gè)步驟:
利用已有的snp數(shù)據(jù)庫，建立相關(guān)性模型券躁，產(chǎn)生重校準(zhǔn)表( recalibration table)
-
根據(jù)這個(gè)模型對原始堿基進(jìn)行調(diào)整惩坑，只會(huì)調(diào)整非已知SNP區(qū)域。
（參考自生信媛--GATK之BaseRecalibrator）
GenomeAnalysisTK -R reference/chr1.fa \
-T BaseRecalibrator \
-knownSites $DATA/dbsnp_138.hg19.vcf \
-o A.recal \
-I A.chr1.sort.dedup.mate.relaign.bam
4）UnifiedGenotyper
這一步提示jimmy的代碼不能用了也拜，說3.7以后采用了新的參數(shù)以舒，更新好快。
使用UnifiedGenotyper尋找variant慢哈。論壇搜索發(fā)現(xiàn)蔓钟，現(xiàn)在推薦用HaplotypeCaller，效果更好卵贱。
GenomeAnalysisTK -T UnifiedGenotyper -R chr1.fa -I A.chr1.sort.dedup.mate.relaign.recal.bam --dbsnp ../common_all_20170710.vcf -o snps.raw.vcf -stand_call_conf 50.0
5）HaplotypeCaller
同樣是更新了的參數(shù)處理方式滥沫。
由于HaplotypeCaller的工作原理，直接省去了BQSR和indel realignment步驟键俱，所以對于一個(gè)variant calling流程而言兰绣，可以直接比對，去重復(fù)后運(yùn)行HaplotypeCaller编振。
HaplotypeCaller缀辩，簡稱HC，能過通過對活躍區(qū)域（也就是與參考基因組不同處較多的區(qū)域）局部重組裝，同時(shí)尋找SNP和INDEL臀玄。換句話說瓢阴，當(dāng)HC看到一個(gè)地方好活躍呀，他就不管之前的比對結(jié)果了健无，直接對這個(gè)地方進(jìn)行重新組裝荣恐，所以就比傳統(tǒng)的基于位置（position-based)的工具，如前任UnifiedGenotyper好用多了累贤。
HC可以處理非二倍體物種和混池實(shí)驗(yàn)數(shù)據(jù)叠穆，但是不推薦用于癌癥或者體細(xì)胞。
HC還可以正確處理可變剪切畦浓，所以也能用于RNA-Seq數(shù)據(jù)痹束。參考自生信媛--生信媛的GATK大禮包
GenomeAnalysisTK -R chr1.fa -T HaplotypeCaller --dbsnp ../common_all_20170710.vcf -stand_call_conf 30 -o A.hc.vcf -I A.chr1.sort.dedup.mate.relaign.recal.bam
4.總結(jié)
雖然馬馬虎虎走了一次最基本的幾個(gè)步驟，發(fā)現(xiàn)學(xué)習(xí)之路任重而道遠(yuǎn)呀讶请！我只是個(gè)初學(xué)者祷嘶，加油！
