二.變異結(jié)果注釋與統(tǒng)計(jì)
書接上回http://www.reibang.com/p/ab6a35502786
如果是使用轉(zhuǎn)錄組重測序的話需要修改過濾條件，會發(fā)現(xiàn)大多數(shù)變異都發(fā)生在基因區(qū)，其他的差別都不大齿梁。比對參考基因組之后進(jìn)行變異檢測姐军，同樣可以用GATK的流程來做。會多一個(gè)cigar字符串的處理步驟汰聋。因?yàn)樵诨蚪M上會比對出跨區(qū)域的reads桑逝，所以要把這樣一整條reads（中間是非基因區(qū)）當(dāng)做兩條來處理弧圆。
三代測序由于錯(cuò)誤率高击费，本身是不適合用于變異檢測的拢蛋，尤其是SNP和InDel。
我們要統(tǒng)計(jì)的SNP變異類型包括：
1.變異位點(diǎn)是否存在于基因區(qū)or上下2kbp范圍內(nèi)蔫巩，如果在2k之內(nèi)的話可能存在一些轉(zhuǎn)錄元件谆棱。太遠(yuǎn)的話就可以視為基因間區(qū)了
2.如果在基因區(qū)的話，就要區(qū)分是否在intron區(qū)圆仔。在UTR區(qū)的變異影響不大垃瞧，而在CDS區(qū)的話又要進(jìn)一步進(jìn)行分析。
3.如果變異影響了氨基酸類型的話坪郭，就會對基因功能（蛋白）造成影響了个从。
4.更加嚴(yán)重的變異類型有：對起始密碼子的突變，突變后形成終止密碼子（提前終止）歪沃，干擾可變剪切位點(diǎn)嗦锐。
在InDel變異中，造成的非3倍數(shù)的移碼突變會非常嚴(yán)重绸罗。即使是3倍數(shù)的移碼突變也是較嚴(yán)重的突變意推。
annovar 注釋
? 文件準(zhǔn)備：
基因組的位置文件gtf
變異位點(diǎn)文件vcf
gtf第八列是相位（phase）信息，意為處在CDS交界處的堿基偏移
SNP 和 INDEL
# 格式轉(zhuǎn)換
$ gtfToGenePred -genePredExt \
genome.gtf \
genome_refGene.txt#下劃線后的東西不能變
# 提取轉(zhuǎn)錄本（來自annovar軟件）
$ retrieve_seq_from_fasta.pl --format refGene \
--seqfile genome.fasta \
--out genome_refGeneMrna.fa \
genome_refGene.txt
輸出得到fa文件之后進(jìn)行注釋珊蟀。annovar可以基于基因位置進(jìn)行注釋，也可以基于指定區(qū)域注釋外驱。
#格式轉(zhuǎn)換
$ convert2annovar.pl \
-format vcf4 \ #輸入文件格式
-allsample \ #不加的話只有第一個(gè)樣本的信息
-withfreq \ #輸出alt頻率
all.filtered.snp.vcf \ #輸入文件育灸，來自GATK的輸出
>all.snp.vcf.annovar.input #輸出文件
#注釋
$ annotate_variation.pl \
-geneanno \ #基于基因進(jìn)行注釋
--neargene 2000 \ #需要考慮的gene上下游范圍
-buildver genome \ #基因組數(shù)據(jù)庫名稱
-dbtype refGene \ #數(shù)據(jù)庫類型
-outfile all.anno.snp \ #輸出文件前綴
-exonsort \ # 結(jié)果按exon進(jìn)行排序
all.snp.vcf.annovar.input ./#最后一項(xiàng)是數(shù)據(jù)庫所在的目錄
#統(tǒng)計(jì)基因組各區(qū)域中snp的個(gè)數(shù)
$ less -N all.anno.variant_function |cut -f 1|sort | uniq -c
分別會得到全基因組變異文件和外顯子變異文件，后者詳細(xì)記載有變異類型昵宇。InDel變異只需要改變文件的名字就可以了磅崭。最后得到的輸出文件中，堿基的突變會不止一個(gè)瓦哎，而變異類型的種類也會更多砸喻。
當(dāng)一個(gè)變異可能屬于多種變異類型的時(shí)候，annovar會按變異造成的后果嚴(yán)重性進(jìn)行優(yōu)先級排序蒋譬，只輸出最嚴(yán)重的割岛。如果需要所有類型的輸出結(jié)果可以用SnpEff軟件或者加參數(shù)。
SV與CNV
可以繼續(xù)用上一步構(gòu)建好的數(shù)據(jù)庫犯助，命令幾乎一樣癣漆。
$ convert2annovar.pl -format vcf4 -allsample -withfreq all_SV.vcf >SV_annovar.input
$ annotate_variation.pl -geneanno --neargene 2000 -buildver genome -dbtype refGene -outfile SV_annovar.output -exonsort SV_annovar.input ./
三.群體結(jié)構(gòu)分析
包括structure，PCA剂买，進(jìn)化樹分析和LD decay等
01.snp數(shù)據(jù)過濾
需要過濾maf次等位基因頻率惠爽，missing癌蓖，還有LD過濾。前二者用于去掉低質(zhì)量snp位點(diǎn)婚肆，后者用于生成符合哈溫平衡的群體租副。
$ plink --vcf all.vcf \ # 輸入文件
--geno 0.1 \ # 設(shè)置缺失率閾值
--maf 0.01 \ # 設(shè)置maf閾值
--biallelic-only strict \ # 去掉三等位位點(diǎn)
--out all.missing_maf \ # 輸出文件前綴
--recode vcf-iid \ # 輸出vcf格式
--allow-extra-chr \ # 允許其他類型的染色體
--set-missing-var-ids @:# \ # 給沒有id的snp取一個(gè)id
--keep-allele-order #讓plink不要給ref和alt換位置
用plink過濾maf和missing，也可以用vcftools來完成较性。
$ plink --vcf all.missing_maf.vcf \
--indep-pairwise 50 10 0.2 \ # LD過濾閾值
--out tmp.ld \#前綴
--allow-extra-chr \
--set-missing-var-ids @:#
$ plink --vcf all.missing_maf.vcf \
--extract tmp.ld.prune.in \ # 保留的SNP列表
--out all.LDfilter \
--recode vcf-iid \#輸出文件的格式
--keep-allele-order \
--allow-extra-chr \
--set-missing-var-ids @:#
第一步輸出是可以保留下來的的位點(diǎn)的id信息附井，所以上一步過濾的時(shí)候最好可以給它加一個(gè)id。第二步再利用位點(diǎn)id信息來提取對應(yīng)的vcf文件两残。
顯然永毅，LD過濾會越來越少。如果不需要那么多snp位點(diǎn)去進(jìn)行分析的話可以把閾值卡得更嚴(yán)格人弓，從而節(jié)省計(jì)算資源沼死。
02.進(jìn)化樹構(gòu)建
? 文件準(zhǔn)備
vcf 文件：all.LDfilter.vcf
$ run_pipeline.pl -Xms1G -Xmx5G \
-SortGenotypeFilePlugin \
-inputFile all.LDfilter.vcf \
-outputFile all.LDfilter.sort.vcf \
-fileType VCF
$ run_pipeline.pl -Xms1G -Xmx5G \
-importGuess all.LDfilter.sort.vcf \
-ExportPlugin \
-saveAs sequences.phy \
-format Phylip_Inter
phylip格式轉(zhuǎn)換之前需要先排序。雖然表面上是一個(gè)pl程序崔赌，但其中有調(diào)用java的部分意蛀，所以仍然需要設(shè)置一些java參數(shù)
dnadist可以用于nj法交互式建樹，不過由于群體遺傳的數(shù)據(jù)量往往非常大健芭，腳本式利用配置文件在后臺運(yùn)行會更加合適县钥。
$ echo -e "sequences.phy\nY" > dnadist.cfg #-e打入換行符
$ dnadist < dnadist.cfg >dnadist.log ##距離矩陣
$ mv outfile infile.dist
$ echo -e "infile.dist\nY" > neighbor.cfg
$ neighbor < neighbor.cfg >nj.log
#距離矩陣
$ less infile.dist | tr '\n' '|'| sed 's/| / /g' | tr '|' '\n' >infile.dist.table
#去掉換行符
$ less outtree | tr '\n' ' '|sed 's/ //g' > outtree.nwk
dnadist發(fā)現(xiàn)outfile會報(bào)錯(cuò)而非覆蓋，所以最好把它名字改掉慈迈。
這樣進(jìn)化樹文件就做好了若贮，再改用其他軟件來作圖。聽說ggtree包很好看痒留。
03.structure分析
# vcf格式轉(zhuǎn)bed格式
$ plink --vcf ../00.filter/all.LDfilter.vcf \
--make-bed --out all --allow-extra-chr \
--keep-allele-order --set-missing-var-ids @:#
# 生成k=2~4的腳本
$ seq 2 4 | \
awk '{print "admixture --cv -j2 all.bed "$1" 1>admix."$1".log 2>&1"}'\
> admixture.sh
$ cat admixture.sh
admixture --cv -j2 all.bed 2 1>admix.2.log 2>&1
admixture --cv -j2 all.bed 3 1>admix.3.log 2>&1
admixture --cv -j2 all.bed 4 1>admix.4.log 2>&1
$ nohup sh admixture.sh &
#繪圖
$ mkdir result
$ cp ./*.Q result/
$ less outtree.nwk |sed 's/:/\n/g'|awk -F '[,(]' '{print $NF}'|awk '$1 !~")"' >sample.pop.order
$ Rscript draw_admixture.R result all.fam sample.pop.order structure
在log文件里面看各K值的CV error
Rscript第一個(gè)參數(shù)是Q矩陣谴麦；第二個(gè)參數(shù)來自格式轉(zhuǎn)換時(shí)生成的fam文件，也可能用別的含樣品順序的文件伸头；第三個(gè)參數(shù)是輸出的樣品順序文件匾效；最后是輸出的前綴。可以基于樹文件中給出的大概分類情況來調(diào)整順序。
其他的看這篇文章吧http://www.reibang.com/p/21c2a683c6fb
04.PCA分析
? 文件準(zhǔn)備：
vcf 文件;all.LDfilter.vcf
分組信息表：sample.pop
plink --vcf all.LDfilter.vcf \
--pca 10 \ # 主成分個(gè)數(shù)
--out PCA_out \
--allow-extra-chr --set-missing-var-ids @:#
利用輸出的PCA_plink.out.eigenvec和分組信息逞姿，在R中進(jìn)行繪圖。如果要分析三個(gè)主成分的話可以輸出兩張圖魔策，一張3dPCA的效果并不好。
05.LDdecay和LDBlock
隨著兩個(gè)site在基因組上的距離變遠(yuǎn)锌妻，連鎖也會變?nèi)醮耍唛g的R2會變小。將所有SNP的LD系數(shù)按距離求平均作為縱軸，將距離作為橫軸
LD衰減距離取決于群體類型搁吓、世代間隔和染色體相對位置.
? 文件準(zhǔn)備
亞群樣品列表文件：
pop.SC.table
pop.YZR.table
$ head -n 3 pop.SC.table
GA0008
GA0035
GA0038
#兩個(gè)亞群
$ PopLDdecay -InVCF all.vcf \ # 輸入vcf文件
-SubPop pop.SC.table \ # 指定要分析的亞群
-MaxDist 500 \ # SNP對距離上限500kb
-OutStat pop.SC.stat
$ PopLDdecay -InVCF ../data/all.vcf -SubPop ../data/pop.YZR.table -MaxDist 500 -OutStat pop.YZR.stat
#多個(gè)亞群共同繪圖
### 準(zhǔn)備配置文件原茅，可以從文件名中awk出來
$ ls pop.*.stat.gz |awk -F"." '{ print $0"\t"$2 }' > ld_stat.list
$ Plot_MultiPop.pl -inList ld_stat.list \
-output ld_stat.multi \
-keepR#用于生成R腳本
LDBlock用于展示染色體局部的相關(guān)性，所以需要自己先有一個(gè)關(guān)注的預(yù)選區(qū)域
$ LDBlockShow -InVCF all.missing_maf.vcf \
-OutPut LDBlockShow \ # 輸出圖片前綴
-Region 1:300000-500000 \ # 繪圖范圍
-SeleVar 2 \ # 繪制R?2
-OutPdf \ # 輸出pdf圖片
-OutPng # 輸出png圖片
LDBlock展示了連鎖遺傳信息堕仔，一般會和曼哈頓圖一起用
四.群體選擇分析
01.π值擂橘，F(xiàn)st，Tajima'D檢驗(yàn)
? 文件準(zhǔn)備
亞群樣品列表文件 pop.SC.table pop.YZR.table 群體變異檢測結(jié)果文件：all.vcf
#π值計(jì)算
$ vcftools --vcf all.vcf \
--window-pi 100000\ #設(shè)置窗口大小
--window-pi-step 10000\ #設(shè)置步長
--keep pop.SC.table \ #亞群樣品列表
--out ./Pi.SC #輸出文件前綴
###同理指定其他亞群
$ vcftools --vcf all.vcf --window-pi 100000 --window-pi-step 10000\ --keep pop.YZR.table --out ./Pi.YZR
#π值沿染色體的分布
$ ls Pi.*pi | awk -F "." '{print $2"\t" $0}' > Pi.list
接下來就可以利用Pi.list來進(jìn)行CMplot繪圖摩骨。染色體中間π值很低的區(qū)域往往是著絲粒區(qū)域通贞。同樣，ROD的曼哈頓圖也可以用CMplot繪制恼五。
#TajimaD
$ vcftools --vcf all.vcf \
--TajimaD 100000 \
--keep ../../data/pop.SC.table \
--out TajimaD.SC#前綴名稱
$ vcftools --vcf all.vcf --TajimaD 100000 --keep pop.YZR.table --out TajimaD.YZR
$ ls TajimaD.*.Tajima.D |awk -F"." '{print $2"\t"$0}' > TajimaD.list
#Fst
$ vcftools --vcf all.vcf \
--fst-window-size $window \ # 窗口大小
--fst-window-step $step \ # 設(shè)置step距離
--weir-fst-pop pop.SC.table \
--weir-fst-pop pop.YZR.table \
--out ./Fst.pop1.pop2
一系列的統(tǒng)計(jì)值計(jì)算都是用vcftools完成的昌罩，在log中有Fst的平均值和加權(quán)值的統(tǒng)計(jì)，可以用加權(quán)平均數(shù)來進(jìn)行CMplot繪圖灾馒。輸出文件為.windowed.weir.fst茎用，兩個(gè)樣本分化越大，這個(gè)值也會越大睬罗。
接下來可以進(jìn)行Fst和ROD的聯(lián)合分析轨功，以及Fst，D和π的聯(lián)合繪圖容达。
02.XP-CLR
這里使用后發(fā)表的Python版本的XPCLR軟件
$ xpclr -I ../data/all.vcf \
-O Chr01.xpclr-python.out \
-C Chr01 \
-Sa pop.SC.table \#兩個(gè)亞群信息
-Sb pop.YZR.table \
--rrate 1e-8 \#單位堿基對應(yīng)的遺傳距離
--ld 0.95 \
--maxsnps 200 \
--size 100000 \#滑窗
--step 20000
$ cut -f 2,3,4,12 Chr01.xpclr-python.out | awk 'NF == 4 ' > Chr01.xpclr-python.table
#提取兩個(gè)chr中xpclr的結(jié)果(第12列）
$ xpclr -I ../data/all.vcf -O Chr02.xpclr-python.out -C Chr02 -Sa pop.SC.table -Sb pop.YZR.table --rrate 1e-8 --ld 0.95 --maxsnps 200 --size 100000 --step 20000
$ cut -f 2,3,4,12 Chr02.xpclr-python.out | awk 'NF == 4 ' > Chr02.xpclr-python.table
#如果是第一行古涧，或者最后一列是xpclr，而且沒有空值的列輸出
$ cat Chr*.xpclr-python.table | awk '{if((NR==1 || $NF != "xpclr") && NF == 4 ){print $0 }}' > all.xpclr-python.table
最后的結(jié)果也可以在R中繪圖花盐，受到選擇越強(qiáng)的區(qū)域xpclr值會越大羡滑。
