一次簡化基因組數(shù)據(jù)分析實(shí)戰(zhàn)

盡管目前已經(jīng)有大量物種基因組釋放出來抒钱，但還是存在許多物種是沒有參考基因組畴蹭。使用基于酶切的二代測序技術(shù)朗恳，如RAD-seq，GBS蔚叨，構(gòu)建遺傳圖譜是研究無參考物種比較常用的方法床蜘。Stacks就是目前比較通用的分析流程，能用來構(gòu)建遺傳圖譜蔑水，處理群體遺傳學(xué)邢锯，構(gòu)建進(jìn)化發(fā)育樹。

這篇教程主要介紹如何使用Stacks分析基于酶切的二代測序結(jié)果搀别，比如說等RAD-seq丹擎，分析步驟為環(huán)境準(zhǔn)備，原始數(shù)據(jù)質(zhì)量評估歇父， 多標(biāo)記數(shù)據(jù)分離蒂培，序列比對（無參則需要進(jìn)行contig de novo 組裝），RAD位點(diǎn)組裝和基因分型榜苫，以及后續(xù)的標(biāo)記過濾和格式轉(zhuǎn)換护戳。

適用范圍：

• 酶切文庫類型：ddRAD, GBS, ezRAD, quad-dRAD和Rapture。 但是stacks更適用于RAD-seq垂睬，GBS推薦TASSEL媳荒。如下是“Genome-wide genetic marker discovery and genotyping using next-generation sequencing”對幾種常見的建庫方法的總結(jié)

幾種文庫的不同

• 測序類型： 雙酶切文庫雙端測序數(shù)據(jù)或單端數(shù)據(jù)
• 測序平臺(tái)： illumina， IonTorren

局限性：

• 不能用于普通的隨機(jī)文庫測序
• 不能適用單酶切文庫的雙端測序數(shù)據(jù)（stacks 2.0可以）
• 無法用于混池測序驹饺，也不適合與多倍體钳枕，因?yàn)閟tacks在組裝時(shí)假定物種為二倍體
• 對于深度不夠，且錯(cuò)誤率比較高的數(shù)據(jù)逻淌，它也沒轍么伯，所以建議深度在20x以上

分析者要求：掌握基本的Unix命令行，會(huì)基本集群操作卡儒，熟悉R語言編程田柔。

硬件要求：電腦的內(nèi)存在64G以上，8～16核CPU起步骨望，準(zhǔn)備1T以上的硬盤硬爆。

前期準(zhǔn)備

準(zhǔn)備分為兩個(gè)部分：軟件安裝和數(shù)據(jù)下載。

數(shù)據(jù)準(zhǔn)備： 數(shù)據(jù)來自于2012年發(fā)表的"The population structure and recent colonization history of Oregon threespine stickleback determined using restriction-site associated DNA-sequencing"中的三刺魚( Gasterosteus aculeatus )數(shù)據(jù)集擎鸠，一共有78個(gè)樣本缀磕，來自于美國俄勒岡海岸的4個(gè)群體，兩個(gè)是海水魚（Cushman Slough’ (CS)和 ‘South Jetty’ (SJ)），兩個(gè)是淡水魚（‘Winchester Creek’ (WC) 和 ‘Pony Creek Reservoir’ (PCR)）袜蚕。

mkdir -p stacks-exercise
cd stacks-exercise
wget -q http://catchenlab.life.illinois.edu/data/rochette2017_gac_or.tar.gz &
tar xf http://catchenlab.life.illinois.edu/data/rochette2017_gac_or.tar.gz

這個(gè)數(shù)據(jù)大約在9G左右糟把，因此需要很長一段時(shí)間，這段時(shí)間可以安裝軟件牲剃。

軟件安裝：需要安裝BWA, SAMtools, stacks,R/ADEgenet. 好消息是這些軟件都可以通過bioconda進(jìn)行安裝遣疯，除了R/ADEgenet推薦直接用R的install.packages("adegenet")

# 適用conda安裝軟件
conda install -c bioconda bwa
conda install -c bioconda samtools
conda install -c bioconda stacks=1.47

估計(jì)此時(shí)數(shù)據(jù)依舊沒有下載完，我們先創(chuàng)建后續(xù)需要用到的文件目錄

mkdir -p stacks-exercise/{00-raw-data,01-clean-data,02-read-alignment,reference/genome,03-stacks-analysis/{de-novo,ref-based},test/{de-novo,ref-based},info}
# 目錄結(jié)構(gòu)
stacks-exercise/
|-- 00-raw-data # 原始數(shù)據(jù)
|-- 01-clean-data # 處理后數(shù)據(jù)
|-- 02-read-alignment # 比對后數(shù)據(jù)
|-- 03-stacks-analysis # 實(shí)際運(yùn)行
|   |-- de-novo
|   |-- ref-based
|-- info # barcode信息
|-- reference # 參考基因組
|   |-- genome
|-- rochette2017_gac_or.tar.gz
|-- test # 測試數(shù)據(jù)集
    |-- de-novo
    |-- ref-based

準(zhǔn)備輸入文件（可選）

這一步并非必須凿傅，取決公司提供給你什么樣的數(shù)據(jù)缠犀。對于多個(gè)樣本測序，公司可能返還的是含有barcode信息原始lane數(shù)據(jù)聪舒，那么就需要從原始數(shù)據(jù)中將各個(gè)樣本的數(shù)據(jù)區(qū)分開辨液。

先將解壓得到的三個(gè)lane的原始數(shù)據(jù)移動(dòng)到我們的文件夾中，

cd stacks-exercise
mv rochette2017_gac_or/top/raw/{lane1,lane2,lane3} 00-raw-data

接著準(zhǔn)備兩個(gè)制表符（Tab）分隔的文件箱残，用于將barcode和樣本對應(yīng)滔迈，以及樣本和群體一一對應(yīng)。這里不需要自己寫了疚宇，只需要將作者存放info里的tsv文件復(fù)制過來即可亡鼠，格式如下

mv rochette2017_gac_or/top/info/*.tsv info/
# barcode和樣本的對應(yīng)關(guān)系
head -n3 info/barcodes.lane1.tsv
CTCGCC  sj_1819.35
GACTCT  sj_1819.31
GAGAGA  sj_1819.32
# 樣本和群體的對應(yīng)關(guān)系
head -n3 info/popmap.tsv
cs_1335.01  cs
cs_1335.02  cs
cs_1335.03  cs

關(guān)barcode和樣本的tsv中，樣本的命名里不要包含空格敷待，只能用字母，數(shù)字仁热，".",“-”和"_", 而且有意義榜揖，最好包含原來群體名的縮寫和樣本編號(hào)。

可視化評估測序數(shù)據(jù)

思考并記錄下按照你的實(shí)驗(yàn)處理抗蠢，你得到的read大概會(huì)是什么結(jié)構(gòu)举哟，從理論上講，從左往右應(yīng)該分別是：barcode迅矛，限制性酶切位點(diǎn)和后續(xù)的測序堿基妨猩。比如說案例應(yīng)該現(xiàn)有6個(gè)堿基的barcode，SbfI限制性位點(diǎn)CCTGCAGG和其余的DNA序列秽褒，總計(jì)101bp

<6-nt barcode>TGCAGG<unique 89-nt sequence>

然后我們就可以用Linux自帶的文本處理命令檢查一下壶硅，比如說grep/zgrep，zcat+head, less/zless.

檢查一下read結(jié)構(gòu)

如果序列已經(jīng)去掉了barcode销斟，那么開頭的就會(huì)是酶切位點(diǎn)庐椒。

第一步：數(shù)據(jù)預(yù)處理

這一步會(huì)用到process_radtags, 它扶負(fù)責(zé)對原始的數(shù)據(jù)進(jìn)行處理，包括樣本分離蚂踊，質(zhì)量控制约谈，檢查酶切位點(diǎn)完整性等。

# 在項(xiàng)目根目錄下
raw_dir=00-raw-data/lane1
barcodes_file=info/barcodes.lane1.tsv
process_radtags -p $raw_dir -b $barcode_file \
    -o 01-clean-data/ -e sbfI --inline_null \
    -c -q -r &> 01-clean-data/process_radtags.lane1.oe &

解釋下參數(shù)，雖然大部分已經(jīng)很明了： -p為原始數(shù)據(jù)存放文件夾棱诱，-b為barcode和樣本對應(yīng)關(guān)系的文件泼橘，-o為輸出文件夾， -e為建庫所用的限制性內(nèi)切酶迈勋，--inline_null表示barcode的位置在單端的read中侥加，-c表示數(shù)據(jù)清洗時(shí)去除表示為N的堿基， -q表示數(shù)據(jù)清理時(shí)要去除低質(zhì)量堿基 -r表示要搶救下barcode和RAD-tag粪躬。

這一步需要留意自己的單端測序担败，還是雙端測序，barcode是在read中镰官，還是在FASTQ的header中提前，是否還需要去接頭序列，是否是雙酶切建庫等泳唠。
另外這一步比較耗時(shí)狈网，盡量脫機(jī)運(yùn)行或者提交到計(jì)算節(jié)點(diǎn)中，不然突然斷網(wǎng)導(dǎo)致運(yùn)行終止就更浪費(fèi)時(shí)間了笨腥。
將運(yùn)行結(jié)果記錄到日志文件中拓哺，方便后期檢查報(bào)錯(cuò)。

運(yùn)行結(jié)束后脖母，在01-clean-data下會(huì)有除了process_radtags.lane1.oe外士鸥，還會(huì)有process_radtags.lane1.log，前者記錄每條lane的數(shù)據(jù)保留情況谆级，后者記錄每個(gè)樣本的數(shù)據(jù)保留情況烤礁。可以將后者復(fù)制到Excel表格中肥照，用柱狀圖等方法直觀了解

樣本剩余read柱狀圖

從圖中可以發(fā)現(xiàn),"sj_1483.05"和"sj_1819.31"幾乎沒有read留下來脚仔，這能是建庫上導(dǎo)致的問題，我們需要將其fastq文件直接刪掉舆绎，從“info/popmap.tsv”中刪掉或者用“#”注釋掉它對應(yīng)行（推薦后者）鲤脏。

在數(shù)據(jù)預(yù)處理這一步，stacks還提供process_shortreads吕朵，clone_filter猎醇， kmer_filter用于處理cDNA文庫和隨機(jī)切割的DNA文庫，如果是RAD-seq不需要用到边锁。

如果是雙端測序姑食，stacks1.47只能通過cat合并兩個(gè)數(shù)據(jù)，而不能有效的利用雙端測序提供的fragment信息茅坛。stacks似乎可以音半，我之后嘗試则拷。

第二步：獲取樣本變異數(shù)據(jù)

這一步之后，分析流程就要根據(jù)是否有參考基因組分別進(jìn)行分析曹鸠。無參考基因組需要先有一步的 de novo 組裝煌茬，產(chǎn)生能用于比對的contig。有參考基因組則需要考慮基因組的質(zhì)量彻桃，如果質(zhì)量太差坛善，則需要進(jìn)一步以無參分析作為補(bǔ)充。

參考基因組主要用于區(qū)分出假陽性的SNP邻眷，將snp與附近其他共線性的snp比較來找出離異值眠屎，這些離異值大多是因?yàn)榻◣爝^程所引入的誤差，如PCR的鏈偏好性擴(kuò)增肆饶。

無論是何者改衩，我們一開始都只能用其中的部分?jǐn)?shù)據(jù)進(jìn)行參數(shù)測試，根據(jù)不同的參數(shù)結(jié)果作為反饋驯镊，進(jìn)行調(diào)優(yōu)葫督，這一步根據(jù)你的運(yùn)氣和經(jīng)驗(yàn)，還有你的算力板惑，時(shí)間不定橄镜。畢竟超算一天，普算一年冯乘。

有參考基因組

三刺魚是可從ensemblgenomic上搜索到到參考基因組信息

• http://asia.ensembl.org/Gasterosteus_aculeatus/Info/Index

但是質(zhì)量非常一般洽胶，僅僅是contig程度，只能說是湊合使用了往湿。

建立索引數(shù)據(jù)庫

stacks不直接參與比對妖异，而是處理不同比對軟件得到的BAM文件，因此你可以根據(jù)自己的喜好選擇比較工具领追。目前，基因組比對工具都是首選BWA-mem响逢，所以這里建立bwa的索引

# 位于項(xiàng)目根目錄下
cd reference/genome
wget -q ftp://ftp.ensembl.org/pub/release-91/fasta/gasterosteus_aculeatus/dna/Gasterosteus_aculeatus.BROADS1.dna.toplevel.fa.gz
gzip -d Gasterosteus_aculeatus.BROADS1.dna.toplevel.fa.gz
cd ..
mkdir -p index/bwa/
genome_fa=genome/Gasterosteus_aculeatus.BROADS1.dna.toplevel.fa
bwa index -p index/bwa/gac $genome_fa &> index/bwa/bwa_index.oe
# 結(jié)果如下
|-- genome
|   |-- Gasterosteus_aculeatus.BROADS1.dna.toplevel.fa
|-- index
    |-- bwa
        |-- bwa_index.oe
        |-- gac.amb
        |-- gac.ann
        |-- gac.bwt
        |-- gac.pac
        |-- gac.sa

小樣本參數(shù)調(diào)優(yōu)

這一步是為了調(diào)整比對工具處理序列相似性的參數(shù)绒窑，保證有絕大多數(shù)的read都能回帖到參考基因組上，因此參數(shù)不能太嚴(yán)格舔亭，能容忍遺傳變異和測序誤差些膨，也不能太寬松，要區(qū)分旁系同源位點(diǎn)钦铺。對于BWA-MEM而言订雾，幾個(gè)和打分相關(guān)的參數(shù)值得注意:

• -B: 不匹配的懲罰, 影響錯(cuò)配數(shù)，默認(rèn)是4
• -O: 缺失和插入的gap打開的懲罰矛洞，影響InDel的數(shù)目洼哎，默認(rèn)是[6,6]
• -E: gap延伸的懲罰烫映，長度k的gap懲罰為'{-O} + {-E}*k'， 默認(rèn)是[1,1]
• -L: soft clip的懲罰噩峦，也就是read兩端直接切掉堿基來保證匹配锭沟，默認(rèn)是[5,5]

對于參考基因組質(zhì)量比較高，且研究物種和參考基因組比較近识补，那么參數(shù)調(diào)整沒有太大必要性族淮。如果質(zhì)量不要，或者所研究物種和參考基因組有點(diǎn)距離凭涂，那么就需要注意不同參數(shù)對結(jié)果的影響祝辣，必要時(shí)使用IGV人工檢查。

讓我們先以默認(rèn)參數(shù)開始切油，處理其中一個(gè)樣本

# 位于項(xiàng)目根目錄下
# 在測試文件下按照參數(shù)創(chuàng)建文件夾
mkdir -p stacks-test/ref-based/{alignment-bwa,stacks-bwa}
## bwa-mem比對
sample=cs_1335.01
fq_file=01-clean-data/$sample.fq.gz
bam_file=test/ref-based/alignment-bwa/${sample}_default.bam
bwa_index=reference/index/bwa/gac
bwa mem -M $bwa_index $fq_file | samtools view -b > $bam_file &

這里的bwa mem使用了-M參數(shù)蝙斜，通常的解釋這個(gè)參數(shù)是為了和后續(xù)的Picard標(biāo)記重復(fù)和GATK找變異兼容。

進(jìn)一步的解釋白翻，不用-M,split read會(huì)被標(biāo)記為SUPPLEMENTARY, 使用該選項(xiàng)則是標(biāo)記為SECONDARY（次要）乍炉，即不是PRIMARY（主要），既然不是主要比對滤馍，所以就被一些工具忽略掉岛琼。如果是SUPPLEMENTARY就有可能在標(biāo)記重復(fù)時(shí)被考慮進(jìn)去。其中split read是嵌合比對的一部分巢株，具體概念見SAM格式解釋槐瑞。

對于比對結(jié)果，可以用samtools stats和samtools flagstat查看一下質(zhì)量

samtools flagstat test/ref-based/alignment-bwa-default/cs_1335.01_default.bam
#1310139 + 0 in total (QC-passed reads + QC-failed reads)
#7972 + 0 secondary
#0 + 0 supplementary
#0 + 0 duplicates
#1271894 + 0 mapped (97.08% : N/A)

97.08%的比對率應(yīng)該是很不錯(cuò)了阁苞，不過可以嘗試降低下錯(cuò)配和gap的懲罰困檩，看看效果

# 位于項(xiàng)目根目錄下
sample=cs_1335.01
fq_file=01-clean-data/$sample.fq.gz
bam_file=test/ref-based/alignment-bwa/${sample}_B3_O55.bam
bwa_index=reference/index/bwa/gac
bwa mem -M -B 3 -O 5,5 $bwa_index $fq_file | samtools view -b > $bam_file &
samtools flagstat test/ref-based/alignment-bwa-default/cs_1335.01_default.bam
#1309830 + 0 in total (QC-passed reads + QC-failed reads)
#7663 + 0 secondary
#0 + 0 supplementary
#0 + 0 duplicates
#1272297 + 0 mapped (97.13% : N/A)

也就提高了0.05%，所以就用默認(rèn)參數(shù)好了那槽。通過IGV可視化悼沿，可以了解簡化基因組的read分布是比較稀疏，10k中可能就只有2個(gè)骚灸。

IGV截圖

得到的BAM文件使用pstack中找snp糟趾，

# 位于項(xiàng)目根目錄下
sample=cs_1335.01
sample_index=1
bam_file=test/ref-based/alignment-bwa/${sample}_B3_O55.bam
log_file=test/ref-based/stacks-bwa/$sample.pstacks.oe
pstacks -t bam -f $bam_file -i $sample_index -o test/ref-based/stacks-bwa/ &> $log_file

這里的參數(shù)也很簡單，-t用來確定輸入文件的格式,-f是輸入文件,-i對樣本編序甚牲，-o指定輸出文件夾义郑。除了以上幾個(gè)參數(shù)外，還可用-p指定線程數(shù)丈钙，-m來制定最低的覆蓋度非驮，默認(rèn)是3.還可以用--model_type [type]制定模型。

最后得到$sample.tags.tsv.gz, $sample.models.tsv.gz, $sample.snps.tsv.gz, 和 $sample.alleles.tsv.gz共4個(gè)文件雏赦，以及一個(gè)日志文件劫笙。參數(shù)評估的主要看日志文件里的幾個(gè)指標(biāo)：

• 實(shí)際使用的alignment數(shù)
• 因soft-clipping剔除的alignment數(shù)芙扎，過高的話要對比對參數(shù)進(jìn)行調(diào)整
• 每個(gè)位點(diǎn)的平均覆蓋度，過低會(huì)影響snp的準(zhǔn)確性邀摆。

這里僅僅用了一個(gè)樣本做測試纵顾，實(shí)際上要用10個(gè)以上樣本做測試，看平均表現(xiàn)栋盹，

全數(shù)據(jù)集處理

在使用小樣本調(diào)試完參數(shù)施逾，這部分參數(shù)就可以應(yīng)用所有的樣本。除了比對和使用pstacks外例获，還需要用到cstacks根據(jù)位置信息進(jìn)一步合并成包含所有位點(diǎn)信息的目錄文件汉额，之后用sstacks從cstacks創(chuàng)建的目錄文件搜索每個(gè)樣本的位點(diǎn)信息。代碼為

cstacks -p 10 --aligned -P 03-stacks-analysis/ref-based -M info/popmap.tsv
# 以其中一個(gè)樣本為例
sample=cs_1335.01
log_file=$sample.sstacks.oe
sstacks --aligned -c 03-stacks-analysis/ref-based/batch_1 -s 03-stacks-analysis/ref-based/$sample -o 03-stacks-analysis/ref-based/ &> 03-stacks-analysis/ref-based/$log_file

你可以寫一個(gè)shell腳本處理榨汤，不過我現(xiàn)在偏好用snakemake寫流程蠕搜，代碼見最后。

無參考基因組

基于參考基因組的分析不能體現(xiàn)出RAD-seq的優(yōu)勢收壕，RAD-seq的優(yōu)勢體現(xiàn)在沒有參考基因組時(shí)他能夠提供大量可靠的分子標(biāo)記妓灌，從而構(gòu)建出遺傳圖譜，既可以用于基因定位蜜宪，也可以輔助組裝虫埂。

和全基因組隨機(jī)文庫不同，RAD-seq是用限制性內(nèi)切酶對基因組的特定位置進(jìn)行切割圃验。這樣的優(yōu)點(diǎn)在于降低了 de novo 組裝的壓力掉伏，原本是根據(jù)overlap（重疊）來延伸150bp左右短讀序列，形成較大的contig澳窑，而現(xiàn)在只是將相似性的序列堆疊(stack)起來斧散。這樣會(huì)產(chǎn)生兩種分子標(biāo)記：1）由于變異導(dǎo)致的酶切位點(diǎn)出現(xiàn)有或無的差異；2）同一個(gè)酶切位點(diǎn)150bp附近存在snp摊聋。

參數(shù)調(diào)優(yōu)

這一步使用的核心工具是ustacks和cstacks,前者根據(jù)序列相似性找出變異鸡捐，后者將變異匯總，同樣需要使用小樣本調(diào)整三個(gè)參數(shù)-M,-n,-m麻裁。 ustacks的 M 控制兩個(gè)不同樣本等位基因（allele）之間的錯(cuò)配數(shù)闯参，m 控制最少需要幾個(gè)相同的堿基來形成一個(gè)堆疊（stack).最后一個(gè)比較復(fù)雜的參數(shù)是能否允許gap(--gap)。而cstacks的 n 和 ustacks的 M等價(jià)悲立。

因此，我們可以嘗試在保證 m=3的前提新博，讓M=n從1到9遞增薪夕，直到找到能讓80%的樣本都有多態(tài)性RAD位點(diǎn)，簡稱r80. Stacks提供了denovo_map.pl來完成這部分工作赫悄，下面開始代碼部分原献。

首先是從info/popmap.tsv中挑選10多個(gè)樣本用于參數(shù)調(diào)試

cat popmap_sample.tsv

cs_1335.01  cs
cs_1335.02  cs
cs_1335.19  cs
pcr_1193.00 pcr
pcr_1193.01 pcr
pcr_1193.02 pcr
pcr_1213.02 pcr
sj_1483.01  sj
sj_1483.02  sj
sj_1483.03  sj
wc_1218.04  wc
wc_1218.05  wc
wc_1218.06  wc
wc_1219.01  wc

然后為每個(gè)參數(shù)都準(zhǔn)備一個(gè)文件夾

mkdir -p test/de-novo/stacks_M{1..9}

然后先 測試 第一個(gè)樣本

# 項(xiàng)目根目錄
M=1
popmap=info/popmap_sample.tsv
reads_dir=01-clean-data
out_dir=test/de-novo/stacks_M${M}
log_file=$out_dir/denovo_map.oe
threads=8
denovo_map.pl -T ${threads} --samples ${reads_dir} -O ${popmap} -o ${out_dir}  -M $M -n $M -m 3 -b 1 -S &> ${log_file} &

代碼運(yùn)行需要一段比較長的時(shí)間馏慨，這個(gè)時(shí)候可以學(xué)習(xí)一下參數(shù)： -T 表示線程數(shù)， --samples表示樣本數(shù)據(jù)所在文件夾姑隅，-O提供需要分析的樣本名写隶， -o是輸出文件夾，之后就是-M, -n, -m這三個(gè)需要調(diào)整的參數(shù)讲仰， -b表示批處理的標(biāo)識(shí)符慕趴， -S關(guān)閉SQL數(shù)據(jù)庫輸出。同時(shí)還有遺傳圖譜和群體分析相關(guān)參數(shù)鄙陡，-p表示為親本冕房，-r表示為后代，-s表明群體各個(gè)樣本趁矾。

為了確保下一步能順利進(jìn)行耙册，還需要對oe結(jié)尾的日志文件的信息進(jìn)行檢查,確保沒有出錯(cuò)

grep -iE "(err|e:|warn|w:|fail|abort)" test/de-novo/stacks_M1/denovo_map.oe

以及以log結(jié)尾的日志中每個(gè)樣本的平均覆蓋度，低于10x的覆蓋度無法保證snp的信息的準(zhǔn)確性

之后對每個(gè)參數(shù)得到的原始數(shù)據(jù)毫捣，要用populations過濾出80%以上都有的snp位點(diǎn)详拙，即r80位點(diǎn)

# 項(xiàng)目根目錄
for M in `seq 1 9`
do
popmap=info/popmap_sample.tsv
stacks_dir=test/de-novo/stacks_M${M}
out_dir=$stacks_dir/populations_r80
mkdir -p ${out_dir}
log_file=$out_dir/populations.oe
populations -P $stacks_dir -O $out_dir -r 0.80 &> $log_file &
done

確定參數(shù)

經(jīng)過漫長的等待后，所有參數(shù)的結(jié)果都已經(jīng)保存到特定文件下蔓同，最后就需要從之確定合適的參數(shù)饶辙，有兩個(gè)指標(biāo)：

• 多態(tài)性位點(diǎn)總數(shù)和r80位點(diǎn)數(shù)
• 短讀堆疊區(qū)的snp平均數(shù)

盡管這些信息都已經(jīng)保存在了相應(yīng)的文本test/de-novo/stacks_M[1-9]/populations_r80/batch_1.populations.log中，但是通過文字信息無法直觀的了解總體情況牌柄，因此更好的方法是用R處理輸出繪圖結(jié)果畸悬。

第一步：提取每個(gè)參數(shù)輸出文件中的log文件中SNPs-per-locus distribution（每個(gè)位點(diǎn)座位SNP分布圖）信息.新建一個(gè)shell腳本，命名為珊佣，log_extractor.sh, 添加如下內(nèi)容

#!/bin/bash
# Extract the SNPs-per-locus distributions (they are reported in the log of populations).
# ----------
echo "Tallying the numbers..."
full_table=n_snps_per_locus.tsv
header='#par_set\tM\tn\tm\tn_snps\tn_loci'
for M in 1 2 3 4 5 6 7 8 9 ;do
    n=$M
    m=3
    # Extract the numbers for this parameter combination.
    log_file=test/de-novo/stacks_M${M}/populations_r80/batch_1.populations.log
    sed -n '/^#n_snps\tn_loci/,/^[^0-9]/ p' $log_file | grep -E '^[0-9]' > $log_file.snps_per_loc
    # Cat the content of this file, prefixing each line with information on this
    # parameter combination.
    line_prefix="M$M-n$n-m$m\t$M\t$n\t$m\t"
    cat $log_file.snps_per_loc | sed -r "s/^/$line_prefix/"
done | sed "1i $header" > $full_table

運(yùn)行后得到n_snps_per_locus.tsv用于后續(xù)作圖蹋宦。會(huì)用到兩個(gè)R腳本plot_n_loci.R和plot_n_snps_per_locus.R，代碼見最后咒锻，結(jié)果圖如下

r80位點(diǎn)數(shù)量

v

SNP分布

從上圖可以發(fā)現(xiàn)M=4以后冷冗，折線就趨于穩(wěn)定，并且每個(gè)座位上的SNP分布趨于穩(wěn)定惑艇，因此選擇M=4作為全樣本數(shù)據(jù)集的運(yùn)行參數(shù)蒿辙。

全數(shù)據(jù)集 de novo 分析

和之前基于參考基因組分析的代碼類型，只不過將序列比對這一塊換成了ustacks滨巴。盡管前面用來確定參數(shù)的腳本denovo_map.pl也能用來批量處理思灌，但它不適合用于大群體分析。ustacks得到的結(jié)果僅能選擇40～200個(gè) 覆蓋度高 且 遺傳多樣性上有代表性的樣本恭取。 使用所有樣本會(huì)帶來計(jì)算上的壓力泰偿，低頻的等位基因位點(diǎn)也并非研究重點(diǎn)，并且會(huì)提高假陽性蜈垮。綜上耗跛，選擇覆蓋度比較高的40個(gè)樣本就行了裕照。

第三步：過濾并導(dǎo)出數(shù)據(jù)

這一步的過濾在stacks1.47是分為兩個(gè)部分。第一部分是對于從頭分析和基于參考基因組都使用rxstacks過濾掉低質(zhì)量變異调塌，然后重新用cstacks和sstacks處理晋南。第二部分是使用population從群體角度進(jìn)行過濾。 在stacks2.0時(shí)代羔砾，rxstacks功能不見了(我推測是作者認(rèn)為作用不大)负间，既然他們都不要了，那我也就只用population過濾和導(dǎo)出數(shù)據(jù)了蜒茄。

這一步主要淘汰那些生物學(xué)上不太合理唉擂，統(tǒng)計(jì)學(xué)上不顯著的位點(diǎn)，一共有四個(gè)主要參數(shù)負(fù)責(zé)控制檀葛，前兩個(gè)是生物學(xué)控制,根據(jù)研究主題而定

• (-r):該位點(diǎn)在單個(gè)群體的所有個(gè)體中的最低比例
• (-p): 該位點(diǎn)至少需要在幾個(gè)群體中存在
• (--min_maf): 過濾過低頻率位點(diǎn)玩祟，推薦5～10%
• (--max_obs_het): 過濾過高雜合率位點(diǎn), 推薦60～70%

# 項(xiàng)目根目錄
min_samples=0.80
min_maf=0.05
max_obs_het=0.70
populations -P 03-stacks-analysis/ref-based/ -r $min_samples --min_maf $min_maf \
--max_obs_het $max_obs_het --genepop &> populations.oe
# --genepop表示輸出格式，可以選擇--vcf等

最后會(huì)在03-stacks-analysis/ref-based/生成至少如下幾個(gè)文件：

• batch_1.sumstats.tsv： 核酸水平上的描述性統(tǒng)計(jì)值屿聋，如觀測值和期望的雜合度 π, and FIS
• batch_1.sumstats_summary.tsv: 每個(gè)群體的均值
• batch_1.hapstats.tsv：單倍型水平的統(tǒng)計(jì)值空扎，如基因多樣性和單倍型多樣性
• batch_1.haplotypes.tsv： 單倍型
• batch_1.genepop：指定的輸出格式
• batch_1.populations.log：輸出日志

至此，上游分析結(jié)束润讥，后續(xù)就是下游分析转锈。后續(xù)更新計(jì)劃：

• stacks總體流程鳥瞰
• Stacks核心參數(shù)深入了解
• RAD-seq和GBS技術(shù)比較
• 不同簡化基因組protocol的比較
• 學(xué)習(xí)TASSEL-GBS數(shù)據(jù)分析流程
• 下游分析探索：這個(gè)得慢慢來

代碼

SAMPLES,       = glob_wildcards("01-clean-data/{sample}.fq.gz")
INDEX_DICT     = {value: key for key, value in dict(enumerate(SAMPLES, start=1)).items()}
FQ_FILES       = expand("01-clean-data/{sample}.fq.gz", sample=SAMPLES)
# de novo
MISMATCH       = 4 # mismatch number of ustacks
DE_NOVO_LOCI   = expand("03-stacks-analysis/de-novo/{sample}.snps.tsv.gz", sample=SAMPLES)
DE_NOVO_CATA   = "03-stacks-analysis/de-novo/batch_1.catalog.snps.tsv.gz"
DE_NOVO_MATS   = expand("03-stacks-analysis/de-novo/{sample}.matches.tsv.gz", sample=SAMPLES)
# ref-based
INDEX          = "reference/index/bwa/gac"
BAM_FILES      = expand("02-read-alignment/ref-based/{sample}.bam", sample=SAMPLES)
REF_BASED_LOCI = expand("03-stacks-analysis/ref-based/{sample}.snps.tsv.gz", sample=SAMPLES)
REF_BASED_CATA = "03-stacks-analysis/ref-based/batch_1.catalog.snps.tsv.gz"
REF_BASED_MATS = expand("03-stacks-analysis/ref-based/{sample}.matches.tsv.gz", sample=SAMPLES)


rule all:
    input: rules.de_novo.input, rules.ref_based.input

rule de_novo:
    input: DE_NOVO_LOCI, DE_NOVO_CATA,DE_NOVO_MATS

rule ref_based:
    input: REF_BASED_LOCI, REF_BASED_CATA, REF_BASED_MATS

# de novo data analysis
## The unique stacks program will take as input a set of short-read sequences
## and align them into exactly-matching stacks (or putative alleles).
rule ustacks:
    input: "01-clean-data/{sample}.fq.gz"
    threads: 8
    params:
        mismatch = MISMATCH,
        outdir   = "03-stacks-analysis/de-novo",
        index    =  lambda wildcards: INDEX_DICT.get(wildcards.sample)
    output:
         "03-stacks-analysis/de-novo/{sample}.snps.tsv.gz",
         "03-stacks-analysis/de-novo/{sample}.tags.tsv.gz",
         "03-stacks-analysis/de-novo/{sample}.models.tsv.gz",
         "03-stacks-analysis/de-novo/{sample}.alleles.tsv.gz",
    log: "03-stacks-analysis/de-novo/{sample}.ustacks.oe"
    shell:"""
    mkdir -p {params.outdir}
    ustacks -p {threads} -M {params.mismatch} -m 3 \
    -f {input} -i {params.index} -o {params.outdir} &> {log}
    """

## choose sample for catalog building
rule choose_representative_samples:
    input: expand("03-stacks-analysis/de-novo/{sample}.ustacks.oe", sample=SAMPLES)
    output:
        "03-stacks-analysis/de-novo/per_sample_mean_coverage.tsv",
        "info/popmap.catalog.tsv"
    shell:"""
    for sample in {input};do
    name=${{sample##*/}}
    name=${{name%%.*}}
    sed -n "/Mean/p" ${{sample}} |\
    sed "s/.*Mean: \(.*\); Std.*/\\1/g" |\
    paste - - - |\
    sed "s/.*/${{name}}\\t&"
    done | sort -k2,2nr > {output[0]}
    head -n 50 {output[0]} | tail -n 40 | cut -f 1 | sed 's/\([0-9a-zA-Z]*\)_.*/&\t\1/' > {output[1]}
    """

rule de_novo_cstacks:
    input:
        "info/popmap.catalog.tsv",
        expand("03-stacks-analysis/de-novo/{sample}.snps.tsv.gz", sample=SAMPLES)
    params:
        stacks_path = "03-stacks-analysis/de-novo/",
        mismatch    = MISMATCH
    threads: 10
    log:"03-stacks-analysis/de-novo/de_novo_cstacks.oe"
    output:
        "03-stacks-analysis/de-novo/batch_1.catalog.alleles.tsv.gz",
        "03-stacks-analysis/de-novo/batch_1.catalog.snps.tsv.gz",
        "03-stacks-analysis/de-novo/batch_1.catalog.tags.tsv.gz"
    shell:"""
        cstacks -p {threads}  -P {params.stacks_path} -M {input[0]} \
        -n {params.mismatch} &> {log}
    """

rule de_novo_sstacks:
    input:
        "03-stacks-analysis/de-novo/{sample}.snps.tsv.gz",
        "03-stacks-analysis/de-novo/{sample}.tags.tsv.gz",
        "03-stacks-analysis/de-novo/{sample}.models.tsv.gz",
        "03-stacks-analysis/de-novo/{sample}.alleles.tsv.gz"
    params:
        catalog = "03-stacks-analysis/de-novo/batch_1",
        sample  = lambda wildcards: "03-stacks-analysis/de-novo/" + wildcards.sample
    output:
        "03-stacks-analysis/de-novo/{sample}.matches.tsv.gz"
    log: "03-stacks-analysis/de-novo/{sample}.sstacks.oe"
    shell:"""
    sstacks  -c {params.catalog} -s {params.sample} -o 03-stacks-analysis/de-novo &> {log}
    """

# reference-based data analysis
## read alignment with bwa-mem
rule bwa_mem:
    input: "01-clean-data/{sample}.fq.gz"
    params:
        index    = INDEX,
        mismatch = "3",
        gap      = "5,5"
    threads: 8
    output: "02-read-alignment/ref-based/{sample}.bam"
    shell:"""
    mkdir -p 02-read-alignment
    bwa mem -t {threads} -M -B {params.mismatch} -O {params.gap} {params.index} {input} \
        | samtools view -b > {output}
    """

## find variant loci
rule pstacks:
    input: "02-read-alignment/ref-based/{sample}.bam"
    params:
        outdir = "03-stacks-analysis/ref-based/",
        index  =  lambda wildcards: INDEX_DICT.get(wildcards.sample)
    threads: 8
    output:
         "03-stacks-analysis/ref-based/{sample}.snps.tsv.gz",
         "03-stacks-analysis/ref-based/{sample}.tags.tsv.gz",
         "03-stacks-analysis/ref-based/{sample}.models.tsv.gz",
         "03-stacks-analysis/ref-based/{sample}.alleles.tsv.gz"
    log: "03-stacks-analysis/ref-based/{sample}.pstacks.oe"
    shell:"""
    mkdir -p 03-stacks-analysis/ref-based
    pstacks -p {threads} -t bam -f {input} -i {params.index} -o {params.outdir} &> {log}
    """

## A catalog can be built from any set of samples processed by the ustacks or pstacks programs
rule ref_based_cstacks:
    input: "info/popmap.tsv",expand("03-stacks-analysis/ref-based/{sample}.snps.tsv.gz", sample=SAMPLES)
    threads: 10
    output:
        "03-stacks-analysis/ref-based/batch_1.catalog.alleles.tsv.gz",
        "03-stacks-analysis/ref-based/batch_1.catalog.snps.tsv.gz",
        "03-stacks-analysis/ref-based/batch_1.catalog.tags.tsv.gz"
    shell:
        "cstacks -p {threads} --aligned -P 03-stacks-analysis/ref-based/ -M {input[0]}"

## Sets of stacks, i.e. putative loci, constructed by the ustacks or pstacks programs
## can be searched against a catalog produced by cstacks.
rule ref_based_sstacks:
    input:
        "03-stacks-analysis/ref-based/{sample}.snps.tsv.gz",
        "03-stacks-analysis/ref-based/{sample}.tags.tsv.gz",
        "03-stacks-analysis/ref-based/{sample}.models.tsv.gz",
        "03-stacks-analysis/ref-based/{sample}.alleles.tsv.gz"
    params:
        catalog = "03-stacks-analysis/ref-based/batch_1",
        sample  = lambda wildcards: "03-stacks-analysis/ref-based/" + wildcards.sample
    output:
        "03-stacks-analysis/ref-based/{sample}.matches.tsv.gz"
    log: "03-stacks-analysis/ref-based/{sample}.sstacks.oe"
    shell:"""
    sstacks --aligned -c {params.catalog} -s {params.sample} -o 03-stacks-analysis/ref-based &> {log}
    """

將以上代碼保存為Snakefile，沒有集群服務(wù)器楚殿，就直接用snakemake運(yùn)行吧撮慨。因?yàn)槲夷茉诩悍?wù)器上提交任務(wù)，所以用如下代碼

snakemake --cluster "qsub -V -cwd" -j 20 --local-cores 10 &

plot_n_snps_per_locus.R的代碼

#!/usr/bin/env Rscript

snps_per_loc = read.delim('./n_snps_per_locus.tsv')
# Keep only M==n, m==3
snps_per_loc = subset(snps_per_loc, M==n & m==3)
# Rename column 1
colnames(snps_per_loc)[1] = 'par_set'

# Create a new data frame to contain the number of loci and polymorphic loci
d = snps_per_loc[,c('par_set', 'M', 'n', 'm')]
d = d[!duplicated(d),]

# Compute these numbers for each parameter set, using the par_set column as an ID
rownames(d) = d$par_set
for(p in rownames(d)) {
    s = subset(snps_per_loc, par_set == p)
    d[p,'n_loci'] = sum(s$n_loci)
    s2 = subset(s, n_snps > 0)
    d[p,'n_loci_poly'] = sum(s2$n_loci)
}

# Make sure the table is ordered
d = d[order(d$M),]

pdf('./n_loci_Mn.pdf')

# Number of loci
# ==========

plot(NULL,
    xlim=range(d$M),
    ylim=range(c(0, d$n_loci)),
    xlab='M==n',
    ylab='Number of loci',
    main='Number of 80%-samples loci as M=n increases',
    xaxt='n',
    las=2
    )
abline(h=0:20*5000, lty='dotted', col='grey50')
axis(1, at=c(1,3,5,7,9))
legend('bottomright', c('All loci', 'Polymorphic loci'), lty=c('solid', 'dashed'))

lines(d$M, d$n_loci)
points(d$M, d$n_loci, cex=0.5)

lines(d$M, d$n_loci_poly, lty='dashed')
points(d$M, d$n_loci_poly, cex=0.5)

# Number of new loci at each step (slope of the previous)
# ==========

# Record the number of new loci at each parameter step
d$new_loci =  d$n_loci - c(NA, d$n_loci)[1:nrow(d)]
d$new_loci_poly = d$n_loci_poly - c(NA, d$n_loci_poly)[1:nrow(d)]

# Record the step size
d$step_size = d$M - c(NA, d$M)[1:(nrow(d))]

plot(NULL,
    xlim=range(d$M),
    ylim=range(c(0, d$new_loci, d$new_loci_poly), na.rm=T),
    xlab='M==n',
    ylab='Number of new loci / step_size (slope)',
    main='Number of new 80%-samples loci as M=n increases'
    )
abline(h=0, lty='dotted', col='grey50')
legend('topright', c('All loci', 'Polymorphic loci'), lty=c('solid', 'dashed'))

lines(d$M, d$new_loci / d$step_size)
points(d$M, d$new_loci / d$step_size, cex=0.5)

lines(d$M, d$new_loci_poly / d$step_size, lty='dashed')
points(d$M, d$new_loci_poly / d$step_size, cex=0.5)

null=dev.off()

plot_n_loci.R的代碼

#!/usr/bin/env Rscript

d = read.delim('./n_snps_per_locus.tsv')
# Keep only M==n, m==3.
d = subset(d, M==n & m==3)
# Make sure the table is ordered by number of snps.
d = d[order(d$n_snps),]

Mn_values = sort(unique(d$M))

# Write the counts in a matrix.
m = matrix(NA, nrow=length(Mn_values), ncol=max(d$n_snps)+1)
for(i in 1:nrow(d)) {
    m[d$M[i],d$n_snps[i]+1] = d$n_loci[i] # [n_snps+1] as column 1 is for loci with 0 SNPs
}

# Truncate the distributions.
max_n_snps = 10
m[,max_n_snps+2] = rowSums(m[,(max_n_snps+2):ncol(m)], na.rm=T)
m = m[,1:(max_n_snps+2)]
m = m / rowSums(m, na.rm=T)

# Draw the barplot.
pdf('n_snps_per_locus.pdf')

clr = rev(heat.colors(length(Mn_values)))

barplot(m,
    beside=T, col=clr, las=1,
    names.arg=c(0:max_n_snps, paste('>', max_n_snps, sep='')),
    xlab='Number of SNPs',
    ylab='Percentage of loci',
    main='Distributions of the number of SNPs per locus\nfor a range of M==n values'
    )
legend('topright', legend=c('M==n', Mn_values), fill=c(NA, clr))

null=dev.off()