三代轉錄組系列：使用Cogent重建基因組編碼區
儘管目前已測序的物種已經很多了，但是對於一些動輒幾個G的複雜基因組，還沒有某個課題組有那麼大的經費去測序，所以仍舊缺少高質量的完整基因組，那麼這個時候一個高質量的轉錄組還是能夠湊合用的。
二代測序的組裝結果只能是差強人意，最好的結果就是不要組裝，直把原模原樣的轉錄組給你是最好的。PacBio Iso-Seq 做的就是這件事情。只不過Iso-Seq測得是轉錄本，由於有些基因存在可變剪下現象，所以所有將一個基因的所有轉錄本放在一起看，搞出一個儘量完整的編碼區。
Cogent能使用高質量的全長轉錄組測序資料對基因組編碼區進行重建，示意圖如下：
軟體安裝
雖然說軟體可以直接通過conda進行安裝，但是根據
下面操作預設你安裝好了miniconda3, 且miniconda3的位置在家目錄下
conda create -n cogent cogent -y
source activate cogent
conda install biopython pulp
conda install networkx==1.10
conda install -c http://conda.anaconda.org/cgat bx-python==0.7.3
pip install -i https://pypi.tuna.tsinghua.edu.cn/simple scikit-image
cd 經過上面一波操作或，請用下面這行命令驗證是否安裝成功
run_mash.py --version
run_mash.py Cogent 3.3
此外還需要安裝另外兩個軟體，分別是Minimap2和Mash
conda install minimap2
wget https://github.com/marbl/Mash/releases/download/v1.0.2/mash-Linux64-v1.0.2.tar.gz
tar xf mash-Linux64-v1.0.2.tar.gz
mv mash $HOME/miniconda3/envs/cogent/bin
對待大資料集，你還需要安裝Cupcake
git clone https://github.com/Magdoll/cDNA_Cupcake.git
cd cDNA_Cupcake
git checkout -b tofu2 tofu2_v21
python setup.py build
python setup.py install
建立偽基因組
讓我們先下載測試所用的資料集，
mkdir test
cd test
https://raw.githubusercontent.com/Magdoll/Cogent/master/test_data/test_human.fa
第一步： 從資料集中搜索同一個基因簇(gene family)的序列
這一步分為超過20,000 條序列的大資料集和低於20,000條序列這兩種情況, 雖然無論那種情況，在這裡我們都只用剛下載的測試資料集
小資料集
第一步：從輸入中計算k-mer譜和配對距離
run_mash.py -k 30 --cpus=12 test_human.fa
# -k, k-mer大小
# --cpus， 程序數
你一定要保證你的輸入是fasta格式，因為該工具目前無法自動判斷輸入是否是fasta格式，所以當你提供詭異的輸入時，它會報錯，然後繼續在後臺折騰。
上面這行命令做的事情是：
- 將輸入的fasta檔案進行拆分，你可以用
--chunk_size指定每個分塊的大小 - 對於每個分塊計算k-mer譜
- 對於每個配對的分塊，計算k-mer相似度（距離）
- 將分塊合併成單個輸出檔案
第二步：處理距離檔案，建立基因簇
process_kmer_to_graph.py test_human.fa test_human.fa.s1000k30.dist test_human human
這裡的test_human是輸出資料夾，human表示輸出檔名字首。此外如果你有輸入檔案中每個轉錄亞型的丰度，那麼你可以用-c引數指定該檔案。
這一步會得到輸出日誌test_human.partition.txt，以及test_human下有每個基因family的次資料夾，檔案裡包含著每個基因簇的相似序列。對於不屬於任何基因family的序列，會在日誌檔案種專門說明，這裡是human.partition.txt
大資料集
如果是超過20,000點大資料集，分析就需要用到Minimap2和Cpucake了。分為如下幾個步驟：
- 使用minimap2對資料進行粗分組
- 對於每個分組，使用上面提到的精細的尋找family工具
- 最後將每個分組的結果進行合併
第一步：使用minimap進行分組
run_preCluster.py要求輸入檔名為isoseq_flnc.fasta, 所以需要先進行軟連線
ln -s test_human.fa isoseq_flnc.fasta
run_preCluster.py --cpus=20
為每個分組構建基因簇尋找命令
generate_batch_cmd_for_Cogent_family_finding.py --cpus=12 --cmd_filename=cmd preCluster.cluster_info.csv preCluster_out test_human
得到的是 cmd 檔案，這個cmd可以直接用bash cmd執行，也可以投遞到任務排程系統。
最後將結果進行合併
printf "Partition\tSize\tMembers\n" > final.partition.txt
ls preCluster_out/*/*partition.txt | xargs -n1 -i sed '1d; $d' {} | cat >> final.partition.txt
第二步：重建編碼基因組
上一步得到每個基因簇, 可以分別重構編碼基因組，所用的命令是reconstruct_contig.py
使用方法: reconstruct_contig.py [-h] [-k KMER_SIZE]
[-p OUTPUT_PREFIX] [-G GENOME_FASTA_MMI]
[-S SPECIES_NAME]
dirname
如果你手頭有一個質量不是很高的基因組，可以使用-G GENOME_FASTA_MMI和-S SPECIES_NAME提供參考基因組的資訊。畢竟有一些內含子是所有轉錄本都缺少的，提供基因組資訊，可以補充這部分缺失。
由於有多個基因簇，所以還需要批量執行命令reconstruct_contig.py
generate_batch_cmd_for_Cogent_reconstruction.py test_human > batch_recont.sh
bash batch_recont.sh
第三步: 建立偽基因組
首先獲取未分配的序列， 這裡用到get_seqs_from_list.py指令碼來自於Cupcake, 你需要將cDNA_Cupcake/sequence新增到的環境變數PATH中。
tail -n 1 human.partition.txt | tr ',' '\n' > unassigned.list
get_seqs_from_list.py hq_isoforms.fasta unassigned.list > unassigned.fasta
這裡的測試資料集裡並沒有未分配的序列，所以上面這一步可省去。
然後將未分配的序列和Cogent contigs進行合併
mkdir collected
cd collected
cat ../test_human/*/cogent2.renamed.fasta ../unassigned.fasta > cogent.fake_genome.fasta
ln -s ../test_human.fa hq_isoforms.fasta
minimap2 -ax splice -t 30 -uf --secondary=no \
cogent.fake_genome.fasta hq_isoforms.fasta > \
hq_isoforms.fasta.sam
sort -k 3,3 -k 4,4n hq_isoforms.fasta.sam > hq_isoforms.fasta.sorted.sam
collapse_isoforms_by_sam.py --input hq_isoforms.fasta -s hq_isoforms.fasta.sorted.sam \
-c 0.95 -i 0.85 --dun-merge-5-shorter \
-o hq_isoforms.fasta.no5merge
由於這裡沒有每個轉錄亞型的丰度檔案cluster_report.csv，所以下面的命令不用執行, 最終結果就是hq_isoforms.fasta.no5merge.collapsed.rep.fa
get_abundance_post_collapse.py hq_isoforms.fasta.no5merge.collapsed cluster_report.csv
filter_away_subset.py hq_isoforms.fasta.no5merge.collapsed
如果運行了上面這行程式碼，那麼最終檔案就應是hq_isoforms.fasta.no5merge.collapsed.filtered.*