1.拿到數(shù)據(jù)后先檢查數(shù)據(jù)是否完整。用md5sum命令杠愧。
#生成md5文件
ls KPGP*| while read KPGP; do echo $KPGP;md5sum ${KPGP} >> ${KPGP}.md5; done
#檢查完整性喳钟，全部顯示OK即可
md5sum -c *.md5
2.對(duì)數(shù)據(jù)進(jìn)行質(zhì)檢。
#質(zhì)檢
nohup fastqc -o /data/XXXX/WGS/01fastqc -t 10 *.fq.gz &
#multiqc合并質(zhì)檢報(bào)告查看
nohup multiqc * -o /data/XXXX/WGS/01fastqc/multiqc &
#若不合格還要進(jìn)行用cutadaptor去接頭等操作减细。此處合格匆瓜，則不贅述。
3.與參考基因組進(jìn)行比對(duì)未蝌。
為什么要比對(duì)驮吱？我們已經(jīng)知道NGS測(cè)序下來的短序列（read）存儲(chǔ)于FASTQ文件里面。雖然它們?cè)径紒碜杂谟行虻幕蚪M萧吠，但在經(jīng)過DNA建庫和測(cè)序之后左冬，文件中不同read之間的前后順序關(guān)系就已經(jīng)全部丟失了。因此纸型，F(xiàn)ASTQ文件中緊挨著的兩條read之間沒有任何位置關(guān)系拇砰，它們都是隨機(jī)來自于原本基因組中某個(gè)位置的短序列而已。
因此狰腌，我們需要先把這一大堆的短序列捋順除破，一個(gè)個(gè)去跟該物種的 參考基因組比較，找到每一條read在參考基因組上的位置癌别，然后按順序排列好皂岔，這個(gè)過程就稱為測(cè)序數(shù)據(jù)的比對(duì)。這 也是核心流程真正意義上的第一步展姐，只有完成了這個(gè)序列比對(duì)我們才有下一步的數(shù)據(jù)分析躁垛。
我們這里將用于流程構(gòu)建的BWA就是其中最優(yōu)秀的一個(gè)剖毯，它將BW(Burrows-Wheeler)壓縮算法和后綴樹相結(jié)合，能夠讓我們以較小的時(shí)間和空間代價(jià)教馆，獲得準(zhǔn)確的序列比對(duì)結(jié)果逊谋。
#bwa建立索引
bwa index -a bwtsw chrom.37.fa
#bwa men
vi bwa.sh
for(( i=1 ; i<=6 ; i++ ))
do
bwa mem -t 4 /data/XXXX/WGS/02hg19/chrom.37.fa \
/data/XXXX/WGS/KPGP-00001_L${i}_R1.fq.gz \
/data/XXXX/WGS/KPGP-00001_L${i}_R2.fq.gz > /data/XXXXX/WGS/03BWA/L${i}.sam
done
nohup sh bwa.sh &
4.sam轉(zhuǎn)為bam，并sort好土铺。
#samtools sam2bam+sort
for(( i=1 ; i<=6 ; i++ ))
do
samtools view L${i}.sam>L${i}.bam
done
for(( i=1 ; i<=6 ; i++ ))
do
samtools sort L${i}.sam L${i}.sort
done
為什么需要排序胶滋？FASTQ文件里面這些被測(cè)序下來的read是隨機(jī)分布于基因組上面的，第一步的比對(duì)是按照FASTQ文件的順序把read逐一定位到參考基因組上之后悲敷，隨即就輸出了究恤，它不會(huì)也不可能在這一步里面能夠自動(dòng)識(shí)別比對(duì)位置的先后位置重排比對(duì)結(jié)果。因此后德，比對(duì)后得到的結(jié)果文件中部宿，每一條記錄之間位置的先后順序是亂的，我們后續(xù)去重復(fù)等步驟都需要在比對(duì)記錄按照順序從小到大排序下來才能進(jìn)行瓢湃。
5.標(biāo)記PCR重復(fù)理张。
為什么要去除這個(gè)duplication？主要是因?yàn)樵赾all snp的時(shí)候绵患，如果某個(gè)變異位點(diǎn)的變異堿基都是來自于PCR重復(fù)雾叭，而我們卻認(rèn)為它深度足夠判斷是真的變異位點(diǎn)，這個(gè)結(jié)論其實(shí)有很大可能是假陽性落蝙。
#安裝好gatk
for(( i=1; i<7; i++ ))
do
gatk MarkDuplicates -I L${i}.sort.bam -O /data/lanyunzhou/WGS/04gatk.mark/L${i}.sort.markdup.bam -M /data/lanyunzhou/WGS/04gatk.mark/L${i}.markdup_metrics.txt
done
#在bwa的時(shí)候沒有添加织狐，因此在這一步添加read group
for(( i=1; i<7; i++ ))
do
java -jar /data/XXXXX/software/software/picard.jar AddOrReplaceReadGroups I=L${i}.sort.markdup.bam O=L${i}.sort.markdup.addrg.bam ID=group1 LB= lib1 PL=illumina PU=unit1 SM=sample1
done
#建立新的bam文件的索引
for(( i=1; i<7; i++ ))
do
samtools index L${i}.sort.markdup.addrg.bam
done
Read Group，是一個(gè)非常重要的信息掘殴，以@RG開頭赚瘦，它是用來將比對(duì)的read進(jìn)行分組的。不同的組之間測(cè)序過程被認(rèn)為是相互獨(dú)立的奏寨。在Read Group中，有如下幾個(gè)信息非常重要：
(1)ID鹰服，這是Read Group的分組ID病瞳，一般設(shè)置為測(cè)序的lane ID（不同lane之間的測(cè)序過程認(rèn)為是獨(dú)立的），下機(jī)數(shù)據(jù)中我們都能看到這個(gè)信息的悲酷，一般都是包含在fastq的文件名中套菜；
(2) PL，指的是所用的測(cè)序平臺(tái)设易，這個(gè)信息不要隨便寫逗柴！特別是當(dāng)我們需要使用GATK進(jìn)行后續(xù)分析的時(shí)候，更是如此顿肺！這是一個(gè)很多新手都容易忽視的一個(gè)地方戏溺，在GATK中渣蜗，PL只允許被設(shè)置為：ILLUMINA，SLX旷祸，SOLEXA耕拷，SOLID，454托享，LS454骚烧，COMPLETE，PACBIO闰围，IONTORRENT赃绊，CAPILLARY，HELICOS或UNKNOWN這幾個(gè)信息羡榴∑敬鳎基本上就是目前市場(chǎng)上存在著的測(cè)序平臺(tái)，
(3) SM炕矮，樣本ID么夫，同樣非常重要，有時(shí)候我們測(cè)序的數(shù)據(jù)比較多的時(shí)候肤视，那么可能會(huì)分成多個(gè)不同的lane分布測(cè)出來档痪，這個(gè)時(shí)候SM名字就是可以用于區(qū)分這些樣本。
(4) LB邢滑，測(cè)序文庫的名字腐螟，這個(gè)重要性稍微低一些，主要也是為了協(xié)助區(qū)分不同的group而存在困后。文庫名字一般可以在下機(jī)的fq文件名中找到乐纸，如果上面的lane ID足夠用于區(qū)分的話，也可以不用設(shè)置LB摇予；
除了以上這四個(gè)之外汽绢，還可以自定義添加其他的信息，不過如無特殊的需要侧戴，對(duì)于序列比對(duì)而言宁昭，這4個(gè)就足夠了。這些信息設(shè)置好之后酗宋，在RG字符串中要用制表符（\t）將它們分開积仗。
bwa mem -t 4 -R ‘@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name’ /path/to/human.fasta read_1.fq.gz read_2.fq.gz > sample_name.sam
6.初步variants calling。（GATK best software for snp+indel）
#生成參考序列索引
samtools faidx /data/XXXX/WGS/02hg19/chrom.37.fa
##建參考序列dict
java -jar /data/XXXXX/software/software/picard.jar CreateSequenceDictionary R=/data/XXXXX/WGS/02hg19/chrom.37.fa O=/data/XXXXX/WGS/02hg19//chrom.37.dict
#用gatk的時(shí)候才可省去局部重比對(duì)蜕猫。用其他軟件不行
#初步calling
for(( i=1; i<7; i++ ))
do
gatk HaplotypeCaller \
-R /data/XXXXX/02hg19/chrom.37.fa \
-I /data/XXXXX/WGS/04gatk.mark/L${i}.sort.markdup.addrg.bam \
-O /data/XXXXX/WGS/05gatk.vcf/L${i}.vcf
done
7.過濾質(zhì)量低的變異寂曹。
for(( i=1; i<6; i++ ))
do
bgzip -f L${i}.vcf
done
for(( i=1; i<6; i++ ))
do
tabix -p vcf L${i}.vcf.gz
done
#SNP
for(( i=1; i<7; i++ ))
do
gatk VariantRecalibrator \
-R /data/XXXXX/WGS/02hg19/chrom.37.fa \
-V /data/XXXXX/WGS/05gatk.vcf/L${i}.vcf.gz \
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 /data/XXXXX/WGS/06VQSR.file/hapmap_3.3.hg19.sites.vcf \
-resource:omini,known=false,training=true,truth=false,prior=12.0 /data/XXXXX/WGS/06VQSR.file/1000G_omni2.5.hg19.sites.vcf \
-resource:1000G,known=false,training=true,truth=false,prior=10.0 /data/XXXXX/WGS/06VQSR.file/1000G_phase1.snps.high_confidence.hg19.sites.vcf \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 /data/XXXXX/WGS/06VQSR.file/dbsnp_138.hg19.vcf \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
-mode SNP \
-tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
--rscript-file /data/XXXXX/WGS/07gatk.snp/L${i}.snps.plots.R \
--tranches-file /data/XXXXX/WGS/07gatk.snp/L${i}.snps.tranches \
-O /data/XXXXX/WGS/07gatk.snp/L${i}.snps.recal
done
for(( i=1; i<7; i++ ))
do
gatk ApplyVQSR \
-R /data/XXXXX/WGS/02hg19/chrom.37.fa \
-V /data/XXXXX/WGS/05gatk.vcf/L${i}.vcf.gz \
-ts-filter-level 99.0 \
--tranches-file /data/XXXXX/WGS/07gatk.snp/L${i}.snps.tranches \
-recal-file /data/XXXXX/WGS/07gatk.snp/L${i}.snps.recal \
-mode SNP \
-O /data/XXXXX/WGS/07gatk.snp/L${i}.snps.VQSR.vcf.gz
done
#INDEL
for(( i=1; i<7; i++ ))
do
gatk VariantRecalibrator \
-R /data/XXXXX/WGS/02hg19/chrom.37.fa \
-O /data/XXXXX/WGS/08gatk.indel/L${i}.snps.indels.recal \
-V /data/XXXXX/WGS/07gatk.snp/L${i}.snps.VQSR.vcf.gz \
-resource:mills,known=true,training=true,truth=true,prior=12.0 /data/XXXXX/WGS/06VQSR.file/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 /data/XXXXX/WGS/06VQSR.file/dbsnp_138.hg19.vcf \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
-mode INDEL \
--max-gaussians 6 \
--rscript-file /data/XXXXX/WGS/08gatk.indel/L${i}.snps.indels.plots.R \
--tranches-file /data/XXXXX/WGS/08gatk.indel/L${i}.snps.indels.tranches
done
for(( i=1; i<7; i++ ))
do
gatk ApplyVQSR \
-R /data/XXXXX/WGS/02hg19/chrom.37.fa \
-V /data/XXXXX/WGS/07gatk.snp/L${i}.snps.VQSR.vcf.gz \
-ts-filter-level 99.0 \
--tranches-file /data/XXXXX/WGS/08gatk.indel/L${i}.snps.indels.tranches \
-recal-file /data/XXXXX/WGS/08gatk.indel/L${i}.snps.indels.recal \
-mode INDEL \
-O /data/XXXXX/WGS/08gatk.indel/L${i}.snps.indel.VQSR.vcf.gz
done
#解壓vcf.gz
bgzip -d *.vcf.gz
8.snpEff對(duì)變異進(jìn)行注釋。
#snpeff (到snpeff.jar那個(gè)文件夾)
java -jar $snpeff GRCh37.75 /data/XXXXX/WGS/07gatk.snp/L1.snps.indel.VQSR.vcf > /data/XXXXX/WGS/09snpeff.note/L1.snp.snpeff.vcf
9.突變頻譜分析（6種）
cat L1.snpeff.vcf |grep -v INDEL |grep -v "^#" |cut -f 1-5 >L1.snps.txt
cat L1.snps.txt |grep -v INDEL |grep -v "^#" |cut -f 4,5|sort |uniq -c |grep -v ",">L1.snps2.txt
dat <- data.frame(type=c('C>A(G>T)','C>T(G>A)','C>G(G>C)','T>A(A>T)','T>G(A>C)','T>C(A>G)'),
counts=c(205713+208287,495660+496650,130847+130034,113079+112752,132902+131814,497236+497094))
library(ggplot2)
p=ggplot(dat,aes(x=type,y=counts,fill=type))+geom_bar(stat="identity")
print(p)
10.祖源分析
祖源分析是用千人基因組數(shù)據(jù)挑位點(diǎn)和自己的基因組數(shù)據(jù)位點(diǎn)進(jìn)行PCA分析。
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4784403/
人類的Y染色體擁有約0.58億個(gè)堿基對(duì)（DNA基本結(jié)構(gòu)）隆圆，約占人類男性體細(xì)胞中DNA的2%漱挚。
人類Y染色體上有86個(gè)基因，這些基因只編碼了23種不同的蛋白質(zhì)匾灶。只有擁有Y染色體才能可能繼承的性狀被稱為雄性性狀棱烂。
人類的Y染色體除了在端粒上的擬常染色體區(qū)的少部分片段（只占有染色體長(zhǎng)度約5%）能與相應(yīng)的X染色體發(fā)生重組，
其外都不能發(fā)生重組阶女。
這些片區(qū)是由原本X染色體與Y染色體同源的片段遺留下來的颊糜。
Y染色體中不能發(fā)生重組的其他部分被稱為“NRY區(qū)”（non-recombining region，非重組區(qū)）秃踩。
這個(gè)區(qū)域中的單核苷酸多態(tài)性被用于父系祖先的追溯衬鱼。
11.SV(CNS)
#delly
delly call -x hg19.excl -o delly.bcf -g hg19.fa input.bam
bcftools view delly.bcf > delly.vcf
- 對(duì)檢測(cè)到的SV進(jìn)行g(shù)enomic feature的注釋
Bioconductor 的intansv包
參考資料：WGS數(shù)據(jù)分析處理流程
