參考基因組沒有,經費也沒那麼多,怎麼辦?

儘管目前已經有大量物種基因組釋放出來，但還是存在許多物種是沒有參考基因組。使用基於酶切的二代測序技術，如RAD-seq，GBS，構建遺傳圖譜是研究無參考物種比較常用的方法。Stacks就是目前比較通用的分析流程，能用來構建遺傳圖譜，處理群體遺傳學，構建進化發育樹。

這篇教程主要介紹如何使用Stacks分析基於酶切的二代測序結果，比如說等RAD-seq，分析步驟為環境準備，原始數據質量評估， 多標記數據分離，序列比對（無參則需要進行contig de novo 組裝），RAD位點組裝和基因分型，以及後續的標記過濾和格式轉換。

適用範圍：

酶切文庫類型：ddRAD, GBS, ezRAD, quad-dRAD和Rapture。 但是stacks更適用於RAD-seq，GBS推薦TASSEL。如下是「Genome-wide genetic marker discovery and genotyping using next-generation sequencing」對幾種常見的建庫方法的總結

幾種文庫的不同

局限性：

不能用於普通的隨機文庫測序

不能適用單酶切文庫的雙端測序數據（stacks 2.0可以）

無法用於混池測序，也不適合與多倍體，因為stacks在組裝時假定物種為二倍體

對於深度不夠，且錯誤率比較高的數據，它也沒轍，所以建議深度在20x以上

分析者要求：掌握基本的Unix命令行，會基本集群操作，熟悉R語言編程。

硬體要求：電腦的內存在64G以上，8～16核CPU起步，準備1T以上的硬碟。

整體分析流程見下圖：

準備分為兩個部分：軟體安裝和數據下載。

數據準備： 數據來自於2012年發表的"The population structure and recent colonization history of Oregon threespine stickleback determined using restriction-site associated DNA-sequencing"中的三刺魚( Gasterosteus aculeatus )數據集，一共有78個樣本，來自於美國俄勒岡海岸的4個群體，兩個是海水魚（Cushman Slough』 (CS)和 『South Jetty』 (SJ)），兩個是淡水魚（『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

這個數據大約在9G左右，因此需要很長一段時間，這段時間可以安裝軟體。

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

# 適用conda安裝軟體

conda install -c bioconda bwa

conda install -c bioconda samtools

conda install -c bioconda stacks=1.47

bioconda怎麼用，

利用conda布署生物信息分析環境

七夕剛過，bioconda中國鏡像來臨

估計此時數據依舊沒有下載完，我們先創建後續需要用到的文件目錄

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}

# 目錄結構

stacks-exercise/

|-- 00-raw-data # 原始數據

|-- 01-clean-data # 處理後數據

|-- 02-read-alignment # 比對後數據

|-- 03-stacks-analysis # 實際運行

| |-- de-novo

| |-- ref-based

|-- info # barcode信息

|-- reference # 參考基因組

| |-- genome

|-- rochette2017_gac_or.tar.gz

|-- test # 測試數據集

 |-- de-novo

 |-- ref-based

這一步並非必須，取決公司提供給你什麼樣的數據。對於多個樣本測序，公司可能返還的是含有barcode信息原始lane數據，那麼就需要從原始數據中將各個樣本的數據區分開。

先將解壓得到的三個lane的原始數據移動到我們的文件夾中，

cd stacks-exercise

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

接著準備兩個制表符（Tab）分隔的文件，用於將barcode和樣本對應，以及樣本和群體一一對應。這裡不需要自己寫了，只需要將作者存放info裡的tsv文件複製過來即可，格式如下

mv rochette2017_gac_or/top/info/*.tsv info/

# barcode和樣本的對應關係

head -n3 info/barcodes.lane1.tsv

CTCGCC sj_1819.35

GACTCT sj_1819.31

GAGAGA sj_1819.32

# 樣本和群體的對應關係

head -n3 info/popmap.tsv

cs_1335.01 cs

cs_1335.02 cs

cs_1335.03 cs

關barcode和樣本的tsv中，樣本的命名裡不要包含空格，只能用字母，數字，".",「-」和"_", 而且有意義，最好包含原來群體名的縮寫和樣本編號。

思考並記錄下按照你的實驗處理，你得到的read大概會是什麼結構，從理論上講，從左往右應該分別是：barcode，限制性酶切位點和後續的測序鹼基。比如說案例應該現有6個鹼基的barcode，SbfI限制性位點CCTGCAGG和其餘的DNA序列，總計101bp

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

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

檢查一下read結構

如果序列已經去掉了barcode，那麼開頭的就會是酶切位點。

這一步會用到 process_radtags, 它扶負責對原始的數據進行處理，包括樣本分離，質量控制，檢查酶切位點完整性等。

# 在項目根目錄下

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 &

解釋下參數，雖然大部分已經很明了： -p為原始數據存放文件夾， -b為barcode和樣本對應關係的文件， -o為輸出文件夾， -e為建庫所用的限制性內切酶， --inline_null表示barcode的位置在單端的read中， -c表示數據清洗時去除表示為N的鹼基， -q表示數據清理時要去除低質量鹼基 -r表示要搶救下barcode和RAD-tag。

這一步需要留意自己的單端測序，還是雙端測序，barcode是在read中，還是在FASTQ的header中，是否還需要去接頭序列，是否是雙酶切建庫等。 另外這一步比較耗時，儘量脫機運行或者提交到計算節點中，不然突然斷網導致運行終止就更浪費時間了。 將運行結果記錄到日誌文件中，方便後期檢查報錯。

運行結束後，在 01-clean-data下會有除了process_radtags.lane1.oe外，還會有process_radtags.lane1.log，前者記錄每條lane的數據保留情況，後者記錄每個樣本的數據保留情況。可以將後者複製到Excel表格中，用柱狀圖等方法直觀了解

樣本剩餘read柱狀圖

從圖中可以發現,"sj_1483.05"和"sj_1819.31"幾乎沒有read留下來，這能是建庫上導致的問題，我們需要將其fastq文件直接刪掉，從「info/popmap.tsv」中刪掉或者用「#」注釋掉它對應行（推薦後者）。

在數據預處理這一步，stacks還提供process_shortreads，clone_filter， kmer_filter用於處理cDNA文庫和隨機切割的DNA文庫，如果是RAD-seq不需要用到。

如果是雙端測序，stacks1.47隻能通過cat合併兩個數據，而不能有效的利用雙端測序提供的fragment信息。stacks似乎可以，我之後嘗試。

這一步之後，分析流程就要根據是否有參考基因組分別進行分析。無參考基因組需要先有一步的 de novo 組裝，產生能用於比對的contig。有參考基因組則需要考慮基因組的質量，如果質量太差，則需要進一步以無參分析作為補充。

參考基因組主要用於區分出假陽性的SNP，將snp與附近其他共線性的snp比較來找出離異值，這些離異值大多是因為建庫過程所引入的誤差，如PCR的鏈偏好性擴增。

無論是何者，我們一開始都只能用其中的部分數據進行參數測試，根據不同的參數結果作為反饋，進行調優，這一步根據你的運氣和經驗，還有你的算力，時間不定。畢竟超算一天，普算一年。

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

但是質量非常一般，僅僅是contig程度，只能說是湊合使用了。

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

# 位於項目根目錄下

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

# 結果如下

|-- genome

| |-- Gasterosteus_aculeatus.BROADS1.dna.toplevel.fa

|-- index

 |-- bwa

 |-- bwa_index.oe

 |-- gac.amb

 |-- gac.ann

 |-- gac.bwt

 |-- gac.pac

 |-- gac.sa

這一步是為了調整比對工具處理序列相似性的參數，保證有絕大多數的read都能回帖到參考基因組上，因此參數不能太嚴格，能容忍遺傳變異和測序誤差，也不能太寬鬆，要區分旁系同源位點。對於BWA-MEM而言，幾個和打分相關的參數值得注意:

-B: 不匹配的懲罰, 影響錯配數，默認是4

-O: 缺失和插入的gap打開的懲罰，影響InDel的數目，默認是[6,6]

-E: gap延伸的懲罰，長度k的gap懲罰為'{-O} + {-E}*k'， 默認是[1,1]

-L: soft clip的懲罰，也就是read兩端直接切掉鹼基來保證匹配，默認是[5,5]

對於參考基因組質量比較高，且研究物種和參考基因組比較近，那麼參數調整沒有太大必要性。如果質量不要，或者所研究物種和參考基因組有點距離，那麼就需要注意不同參數對結果的影響，必要時使用IGV人工檢查。

讓我們先以默認參數開始，處理其中一個樣本

# 位於項目根目錄下

# 在測試文件下按照參數創建文件夾

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參數，通常的解釋這個參數是為了和後續的Picard標記重複和GATK找變異兼容。

進一步的解釋，不用 -M,split read會被標記為SUPPLEMENTARY, 使用該選項則是標記為SECONDARY（次要），即不是PRIMARY（主要），既然不是主要比對，所以就被一些工具忽略掉。如果是SUPPLEMENTARY就有可能在標記重複時被考慮進去。其中split read是嵌合比對的一部分，具體概念見SAM格式解釋。

對於比對結果，可以用 samtools stats和 samtools flagstat查看一下質量

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%的比對率應該是很不錯了，不過可以嘗試降低下錯配和gap的懲罰，看看效果

# 位於項目根目錄下

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%，所以就用默認參數好了。通過IGV可視化，可以了解簡化基因組的read分布是比較稀疏，10k中可能就只有2個。

IGV截圖

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

# 位於項目根目錄下

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

這裡的參數也很簡單， -t用來確定輸入文件的格式, -f是輸入文件, -i對樣本編序， -o指定輸出文件夾。除了以上幾個參數外，還可用 -p指定線程數， -m來制定最低的覆蓋度，默認是3.還可以用 --model_type [type]制定模型。

最後得到 $sample.tags.tsv.gz, $sample.models.tsv.gz, $sample.snps.tsv.gz, 和 $sample.alleles.tsv.gz共4個文件，以及一個日誌文件。參數評估的主要看日誌文件裡的幾個指標：

這裡僅僅用了一個樣本做測試，實際上要用10個以上樣本做測試，看平均表現，

全數據集處理在使用小樣本調試完參數，這部分參數就可以應用所有的樣本。除了比對和使用pstacks外，還需要用到 cstacks根據位置信息進一步合併成包含所有位點信息的目錄文件，之後用 sstacks從 cstacks創建的目錄文件搜索每個樣本的位點信息。代碼為

cstacks -p 10 --aligned -P 03-stacks-analysis/ref-based -M info/popmap.tsv

# 以其中一個樣本為例

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

你可以寫一個shell腳本處理，不過我現在偏好用snakemake寫流程，代碼見最後。

基於參考基因組的分析不能體現出RAD-seq的優勢，RAD-seq的優勢體現在沒有參考基因組時他能夠提供大量可靠的分子標記，從而構建出遺傳圖譜，既可以用於基因定位，也可以輔助組裝。

和全基因組隨機文庫不同，RAD-seq是用限制性內切酶對基因組的特定位置進行切割。這樣的優點在於降低了 de novo 組裝的壓力，原本是根據overlap（重疊）來延伸150bp左右短讀序列，形成較大的contig，而現在只是將相似性的序列堆疊(stack)起來。這樣會產生兩種分子標記：1）由於變異導致的酶切位點出現有或無的差異；2）同一個酶切位點150bp附近存在snp。

這一步使用的核心工具是 ustacks和 cstacks,前者根據序列相似性找出變異，後者將變異匯總，同樣需要使用小樣本調整三個參數 -M, -n, -m。 ustacks的 M 控制兩個不同樣本等位基因（allele）之間的錯配數，m 控制最少需要幾個相同的鹼基來形成一個堆疊（stack).最後一個比較複雜的參數是能否允許gap(--gap)。而 cstacks的 n 和 ustacks的 M等價。

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

首先是從 info/popmap.tsv中挑選10多個樣本用於參數調試

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

然後為每個參數都準備一個文件夾

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

然後先 測試 第一個樣本

# 項目根目錄

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} &

代碼運行需要一段比較長的時間，這個時候可以學習一下參數： -T 表示線程數， --samples表示樣本數據所在文件夾，-O提供需要分析的樣本名， -o是輸出文件夾，之後就是-M, -n, -m這三個需要調整的參數， -b表示批處理的標識符， -S關閉SQL資料庫輸出。同時還有遺傳圖譜和群體分析相關參數，-p表示為親本，-r表示為後代，-s表明群體各個樣本。

為了確保下一步能順利進行，還需要對oe結尾的日誌文件的信息進行檢查,確保沒有出錯

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

以及以log結尾的日誌中每個樣本的平均覆蓋度，低於10x的覆蓋度無法保證snp的信息的準確性

之後對每個參數得到的原始數據，要用 populations過濾出80%以上都有的snp位點，即r80位點

# 項目根目錄

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

經過漫長的等待後，所有參數的結果都已經保存到特定文件下，最後就需要從之確定合適的參數，有兩個指標：

多態性位點總數和r80位點數

短讀堆疊區的snp平均數

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

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

#!/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

運行後得到 n_snps_per_locus.tsv用於後續作圖。會用到兩個R腳本 plot_n_loci.R和 plot_n_snps_per_locus.R，代碼見最後，結果圖如下

r80位點數量SNP分布

從上圖可以發現M=4以後，折線就趨於穩定，並且每個座位上的SNP分布趨於穩定，因此選擇M=4作為全樣本數據集的運行參數。

和之前基於參考基因組分析的代碼類型，只不過將序列比對這一塊換成了 ustacks。儘管前面用來確定參數的腳本 denovo_map.pl也能用來批量處理，但它不適合用於大群體分析。 ustacks得到的結果僅能選擇40～200個 覆蓋度高 且 遺傳多樣性上有代表性的樣本。 使用所有樣本會帶來計算上的壓力，低頻的等位基因位點也並非研究重點，並且會提高假陽性。綜上，選擇覆蓋度比較高的40個樣本就行了。

這一步的過濾在stacks1.47是分為兩個部分。第一部分是對於從頭分析和基於參考基因組都使用 rxstacks過濾掉低質量變異，然後重新用 cstacks和 sstacks處理。第二部分是使用 population從群體角度進行過濾。 在stacks2.0時代， rxstacks功能不見了(我推測是作者認為作用不大)，既然他們都不要了，那我也就只用 population過濾和導出數據了。

這一步主要淘汰那些生物學上不太合理，統計學上不顯著的位點，一共有四個主要參數負責控制，前兩個是生物學控制,根據研究主題而定

# 項目根目錄

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等

最後會在 03-stacks-analysis/ref-based/生成至少如下幾個文件：

batch_1.sumstats.tsv： 核酸水平上的描述性統計值，如觀測值和期望的雜合度 π, and FIS

batch_1.sumstats_summary.tsv: 每個群體的均值

batch_1.hapstats.tsv：單倍型水平的統計值，如基因多樣性和單倍型多樣性

batch_1.haplotypes.tsv： 單倍型

batch_1.genepop：指定的輸出格式

batch_1.populations.log：輸出日誌

至此，上遊分析結束，後續就是下遊分析。後續更新計劃：

stacks總體流程鳥瞰

Stacks核心參數深入了解

RAD-seq和GBS技術比較

不同簡化基因組protocol的比較

學習TASSEL-GBS數據分析流程

下遊分析探索：這個得慢慢來

snakfile

關於snakemake的用法，vip論壇有比較好的筆記，春節中獎的小夥伴記得多逛逛

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，沒有集群伺服器，就直接用 snakemake運行吧。因為我能在集群伺服器上提交任務，所以用如下代碼

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()