將經(jīng)過(guò)過(guò)濾后得到的高質(zhì)量數(shù)據(jù)比對(duì)到參考基因組上,并進(jìn)行排序和去重復(fù)等處理,用于后續(xù)的變異檢測(cè)。
軟件準(zhǔn)備
bwa
samtools
picard
文件準(zhǔn)備
參考基因組:genome.fasta
測(cè)序數(shù)據(jù)文件:sample1_1.fq.gz、sample1_2.fq.gz、sample2_1.fq.gz、sample2_2.fq.gz
構(gòu)建基因組index
重測(cè)序很多軟件不能直接讀取基因組,需要對(duì)基因組構(gòu)建index才能使用。
samtools faidx genome.fasta #基因組的統(tǒng)計(jì)信息
bwa index genome.fasta #bwa比對(duì)需要的索引文件
picard CreateSequenceDictionary R=genome.fasta #生成字典序文件
比對(duì)
采用bwa (Li,H.,and R. Durbin等,2009) mem程序?qū)⒔?jīng)過(guò)過(guò)濾后得到的高質(zhì)量數(shù)據(jù)比對(duì)到參考基因組上,比對(duì)的參數(shù)均按照bwamem的默認(rèn)參數(shù)。采用samtools軟件對(duì)sam文件進(jìn)行排序并轉(zhuǎn)換為bam文件。測(cè)序數(shù)據(jù)的生成過(guò)程涉及到文庫(kù)的擴(kuò)增和簇的形成,這兩個(gè)步驟容易產(chǎn)生一些 Duplicates(即PCR Duplicates和Optical Duplicates),這些Duplicates 不能作為變異檢測(cè)的證據(jù)。采用Picard軟件包中的“MarkDuplicates” 去除Duplicates,即如果多個(gè)Paired Reads比對(duì)后具有相同的染色體坐標(biāo),則僅保留分值最高的Paired Reads。
bwa mem -t 2 \ # -t 線(xiàn)程數(shù)
-R '@RG\tID:{sample}\tSM:{sample}\tPL:illumina' \ # 添加表頭信息
genome.fasta \ 參考基因組index路徑
S1_1.fq.gz \ fq1數(shù)據(jù)
S1_2.fq.gz \ fq2數(shù)據(jù)
2>S1.bwa.log | \ 寫(xiě)到日志文件
samtools sort -@ 2 -m 1G -o S1.sort.bam - \ 排序
##生成文件:S1.sort.bam
#去除重復(fù)
picard -Xmx4g \ #設(shè)置內(nèi)存
MarkDuplicates I=S1.sort.bam \ #輸入文件
O=S1.sort.rmdup.bam \ #輸出文件
CREATE_INDEX=true \ #創(chuàng)建bam index
REMOVE_DUPLICATES=true M=S1.marked_dup_metrics.txt
#比對(duì)率統(tǒng)計(jì)
samtools flagstat S1.sort.bam > S1.sort.bam.flagstat
#覆蓋率統(tǒng)計(jì)
samtools coverage S1.sort.bam