前前後後接觸了一些基因組和轉錄組拼接的工作，而且後期還會持續進行。期間遇到了各種各樣莫名其妙的坑，也嘗試了一些不同的方法和軟體，簡單做一個階段性小結。上週的今天更新了原理部分 二代測序資料拼接之原理篇 (點選閱讀)，本篇是閹割版應用部分（原文程式碼太多影響閱讀體驗）。

10000 字(含程式碼)，約20分鐘，文|思考問題的熊

感覺本文過長可直接收藏，然後閱讀 2018，從“丟”開始（點選閱讀）

拼接大致流程

流程的前面4步和DBG演算法相關，尤其是一二兩步需要較大的記憶體。拼接結果受 kmer size，kmer coverage cutoff 和 length and coverage parameters 的影響

資料預處理

去接頭和低質量reads

類似於通常 RNA-seq 資料處理，略~

使用軟體 khmer 進行標準化

digital normalization

關於是否進行 digital normalization 其實一直有不少討論，最初一篇介紹 digital normalization 的文章指出，所謂digital normalization，是 discards redundant data and both sampling variation and the number of errors present in deep sequencing data sets。其最大的好處是可以降低拼接對記憶體的要求並且節省時間，而且對於拼出的 contig 沒有什麼影響。之所以不影響拼接質量，是因為並沒有去掉那些低覆蓋度的資料。具體可以參考這篇文章 What is digital normalization, anyway。

當然，過了一段時間，還是上文作者又寫了一篇部落格，說明digital normalization 存在的一些問題。

另附 軟體官網使用說明

軟體主要功能

• normalizing read coverage ("digital normalization")
• dividing reads into disjoint sets that do not - connect ("partitioning")
• eliminating reads that will not be used by a de - Bruijn graph assembler;
• removing reads with low- or high-abundance k-mers;
• trimming reads of certain kinds of sequencing errors;
• counting k-mers and estimating data set coverage based on k-mer counts;
• running Velvet and calculating assembly statistics;
• converting FASTQ to FASTA;
• converting between paired and interleaved formats for paired FASTQ data

在使用khmer處理資料之前，首先要想清楚是否進行處理。其中最重要的引數是 graph-size filtering 和 graph partitioning。這個軟體拼接的時候可以用，計算表達量差異的切忌使用。

另一個是關於 memory 的設定問題，在官方給出的建議中說了一大堆，總的來說就是越大越好 :)

建議使用伺服器總記憶體的一半，如果記憶體不夠的話會報錯。一般1 billion 的 mRNAseq 需要 -M16G 16G記憶體。如果kmer過小，在進行資料清洗的時候，很可能會造成誤傷。

Khmer 的四種用法

• k-mer counting and abundance filtering
• Partitioning
• Digital normalization
• Read handling: interleaving, splitting, etc.

這裡主要使用Digital normalization

關於kmer設定的說明

The interaction between these three parameters and the filtering process is complex and depends on the data set being processed, but higher coverage levels and longer k-mer sizes result in less data being removed. Lower memory allocation increases the rate at which reads are removed due to erroneous estimates of their abundance, but this process is very robust in practice

針對mRANseq的拼接，官方文獻給出的建議是

By normalizing to a higher coverage of 20, removing errors, and only then reducing coverage to 5, digital normalization can retain accurate reads for most regions. Note that this three-pass protocol is not considerably more computationally expensive than the single-pass protocol: the first pass discards the majority of data and errors, so later passes are less time and memory intensive than the first pass

但是針對本身數量不是過大的資料，目前推薦使用single-pass digital normalization的方法。

pair end 資料合併

for i in `ls /projects/zhaofei/wheat_rawdata/LF0*_1.fq.gz`
do
id=`basename $i |sed 's/_1.fq.gz//'`
interleave-reads.py $i /projects/zhaofei/wheat_rawdata/${id}_2.fq.gz -o ${id}_pair.fq.gz --gzip
done

for i in `ls /projects/zhaofei/wheat_rawdata/LF1*_1.fq.gz`
do
id=`basename $i |sed 's/_1.fq.gz//'`
interleave-reads.py $i /projects/zhaofei/wheat_rawdata/${id}_2.fq.gz -o ${id}_pair.fq.gz --gzip
done

single-pass digital normalization

核心步驟，設定相應的cutoff和kmer進行資料處理，cutoff 20； kmer 21/25/27/31/

此處的cutoff 是指：when the median k-mer coverage level is above this number the read is not kept

normalize-by-median.py -p -k 27 -M 50G -C 20 -R LF01.log 
-o LF01k32c28.fq.gz --gzip LF01_pair.fq.gz > runLF01.log 2>&1 &

如果在去街頭和低質量資料過程中產生了單端資料，可以在命令中加入-u se.fq 引數

如果想保留生成的dbg檔案，可以新增引數 --savegraph normC20k27.ct

去除可能錯誤的kmer

這一步去除上一步中coverage很高，但是kmer abundanc 低的reads。

filter-abund.py -V -Z 18 normC20k20.ct input.keep.fa && 
  rm input.keep.fa normC20k20.ct

# 我自己使用的時候實際沒有執行這一步

提取pair end reads

extract-paired-reads.py input.keep.fq

這一步會分別生成仍是pair reads和非 pair reads，生成的資料可以用來後續正式的拼接過程。

分離pair end reads

split-paired-reads.py input.fq.pe

分離後得到的兩個fastq檔案就可以正式的拼接了。

使用 Trinity 進行拼接

軟體介紹

下載地址：

https://github.com/trinityrnaseq/trinityrnaseq/wiki

相關文獻：

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3571712/

Trinity是目前最常用的轉錄組拼接軟體。

拼接過程共分為三步

Inchworm: 拼接過程。從reads到contigs的過程，contigs可能代表一個全長轉錄本或者一個轉錄本的一部分。會提取所有的重疊k-mers，根據丰度高低檢查每個k-mers，然後將重疊的k-mers延長，形成一個contig

Chrysalis: 將上一步生成的contigs聚類，對每個聚類結果構建DBG圖。一個DBG圖代表了一個基因的全長轉錄本。

Butterfly: 根據上一步構建的DBG圖和圖中的pair end reads 資訊尋找最優路徑。得到具有可變剪接的全長轉錄本，同時分開旁系基因的轉錄本。

直系同源的序列因物種形成（speciation）而被區分開（separated）：若一個基因原先存在於某個物種，而該物種分化為了兩個物種，那麼新物種中的基因是直系同源的；旁系同源的序列因基因複製（gene duplication）而被區分開（separated）：若生物體中的某個基因被複制了，那麼兩個副本序列就是旁系同源的。

常用命令示例

trinityrnaseq-Trinity-v2.4.0/Trinity 
--seqType fq --max_memory 40G --CPU 4 
--left left.fq.gz 
--right right.fq.gz 
--output test_trinity --full_cleanup 
--no_version_check > test.log 2>&1 &

想要詳細瞭解trinity的處理過程，只需要認真讀一下生成的log檔案就可以。

大體上是如下幾步：

---------------------------
Trinity Phase 1: Clustering of RNA-Seq Reads
---------------------------

In silico Read Normalization
-- (Removing Excess Reads Beyond 50 Coverage --

Jellyfish
-- (building a k-mer catalog from reads) --

Inchworm
-- (Linear contig construction from k-mers) --

Chrysalis
-- (Contig Clustering & de Bruijn Graph Construction) --

------------------------
 Trinity Phase 2: Assembling Clusters of Reads
 ---------------------

Butterfly assemblies are written to 
/projects/zhaofei/wheat_assembly/trinity/LF20_1_trinity.Trinity.fasta

可能出現的報錯

需要注意的是，有時候使用trinity拼接一些公用資料會報錯

#If your data come from SRA, be sure to dump the fastq file like so:
 #SRA_TOOLKIT/fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files --gzip file.sra

可以使用命令

fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files --gzip input.sra

有的時候不是公用資料仍然報錯，是因為pair end 資料第一行@開頭的名字不能有空格，必須用1/;2/結尾

#如果中間有空格，結尾正確

 zcat ERR037683_1.fastq.gz|awk '{ if (NR%4==1) { print $1"_"$2"" } else { print } }'|gzip >ERR037683_new_1.fastq.gz

 #跑迴圈修改
 for i in ERR037679_2.fastq.gz ERR037681_2.fastq.gz ERR037687_2.fastq.gz
  do
  zcat $i | awk '{ if (NR%4==1) { print $1"_"$2"" } else { print } }'|gzip > ${i/_2.fastq.gz/}_new_2.fastq.gz && rm -f $i
 done

 #如果結尾不正確
 zcat ERR037683_1.fastq.gz|awk '{ if (NR%4==1) { print $1"_"$2"/1" } else { print } }'|gzip >ERR037683_new_1.fastq.gz
 #_2.fq 同理替換

測試小眾Bridger拼接軟體

github 地址：

https://github.com/fmaguire/Bridger_Assembler

文獻

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0596-2

按照文章裡的這個軟體結合了trinity和cufflinks的優點，和trinity相比拼接速度更快佔記憶體更小，可以產生更好的contig(注意：不一定是好事)

安裝軟體

首先安裝boost

a) download latest boost and unpack it.

        $ tar zxvf boost_1_47_0.tar.gz

     b) change to the boost directory and run ./bootstrap.sh.

        $ cd  boost_1_47_0
        $ ./bootstrap.sh

        $ ./b2 install --prefix=<YOUR_BOOST_INSTALL_DIRECTORY>

        For example,if you want install boost in /home/czheng/local/boost,the commnd is :
        $ ./b2 install --prefix=/home/czheng/local/boost

        If the boost is installed successfully, you would fild two sub-directories in /home/czheng/local/boost/:
        /home/czheng/local/boost/include/
        /home/czheng/local/boost/lib/

       Note: The default Boost installation directory is /usr/local. Take note of the boost installation
        directory, beacuase you need to tell the Bridger installer where to find boost later on.

     d) Set the LD_LIBRARY_PATH enviroment variable:

        The ~/.bash_profile ($HOME/.bash_profile) or ~/.profile file is executed when you login using console or remotely using ssh.
        Append the following command to  ~/.bash_profile or ~/.profile file:
        $ export LD_LIBRARY_PATH=/home/czheng/local/boost/lib:$LD_LIBRARY_PATH

        Save and close the file.

        OR

        just type the command:
        $ export LD_LIBRARY_PATH=/home/czheng/local/boost/lib:$LD_LIBRARY_PATH

        Note: please replace "/home/czheng/local/boost/lib" with your own directory "<YOUR_BOOST_INSTALL_DIRECTORY>/lib"
        If you do not set this variable , you would possible see the follwoing error information:
         "error while loading shared libraries: libboost_serialization.so.1.47.0: cannot open shared object file: No such file or dire    ctory"

再安裝bridger

Building Bridger [Make sure Boost has been installed successfully]

     a) Unpack the Bridger and change to the Bridger direcotry.

        $ tar zxvf Bridger_r2013-06-02.tar.gz
        $ cd Bridger_r2013-06-02

     b) Configure Bridger. If Boost is installed somewhere other than /usr/local, you will need to tell
        the installer where to find it using --with-boost option.

        $ ./configure --with-boost=/home/czheng/local/boost/
        Note: please replace "/home/czheng/local/boost/" with your own directory "<YOUR_BOOST_INSTALL_DIRECTORY>"


     c) Make Bridger.

        $ make

       note: If you build boost suffessfully without using --prefix option, the following commands may need before "make":
           export LIBS="-L/home/czheng/boost_1_47_0/stage/lib" (replace "/home/czheng/boost_1_47_0/" with your own directory)
           export CPPFLAGS="-I/home/czheng/boost_1_47_0/"

報錯和解決方法

如果安裝最新版本的軟體，是會報錯的。這個問題已經有人在GitHub上提出來了，但是軟體的作者已經到阿里巴巴上班了，於是只能自己想辦法改一下原始程式碼。

Splicing Graphs Reconstruction
CMD: ~/software/Bridger_r2014-12-01/src/Assemble --reads both.fa -k 25 --pair_end --fr_strand 2 2>Assemble.log
Error, cmd: ~/software/Bridger_r2014-12-01/src/Assemble --reads both.fa -k 25 --pair_end --fr_strand 2 2>Assemble.log died with ret 256 !

Assemble.log:
[Error] Cannot create directory ./RawGraphs/ !

出現報錯以後，如果檢視目錄會發現已經有/RawGraphs/目錄，嘗試更改相關原始碼解決問題。比如增加生成資料夾的一步。

修改 Bridger.pl 檔案

## Assemble step:
print "n### Splicing Graphs Reconstruction ###nn";
my $graphdir = "RawGraphs"; #新加
mkdir "$output_directory/$graphdir"; #新加
my $assemble_cmd = "$SRC_DIR/Assemble --reads $target_fa -k $kmer_length ";

核心命令

perl Bridger_r2014-12-01/Bridger.pl --seqType fq -k 31 --left left.fq --right right.fq --CPU 6 --debug -o output

cd-hit 聚類去冗餘

cd-hit 是一個聚類軟體，可以對DNA序列或者蛋白質序列進行聚類，其本質應該還是多序列比對。為了讓拼出來的scaffold更少一些，可以嘗試對拼接結果再次聚類，輸出每個聚類結果中最長的序列。類似的軟體還有corset和lace。

但是實際使用過程中，發現效果也一定會特別明顯。

cd-hit-est -i input.fasta -o output-cdhit -T 10 -M 200000

會輸出兩個結果檔案，一個包含序列資訊，一個是fasta檔案。

評估拼接質量

使用軟體 transrate

相關文獻 ：http://genome.cshlp.org/content/early/2016/06/01/gr.196469.115

官方網站 ：http://hibberdlab.com/transrate/index.html

主要有三個功能

• by inspecting the contig sequences
• by mapping reads to the contigs and inspecting the alignments
• by aligning the contigs against proteins or transcripts from a related species and inspecting the alignments

這裡我們主要使用前兩個功能，如果是有參轉錄組的拼接，可以嘗試使用第三個。但如果是為了檢視新的轉錄本，進行第三項評估也沒有太大意義。針對轉錄組拼接而言，第一步中各種長度的統計結果意義也不大，只有回帖率這個指標是最重要的。通過第二部評估，transrate會返回非常多的有用資訊。具體結果解讀可以參考網站。需要說明的是，在輸出結果中，會直接生成一個good.fasta 檔案，本質是在計算轉錄本表達量時不為0的序列。

這個軟體進行第二部計算時，呼叫了SNAP 進行map，呼叫了salmon 評估轉錄本表達量，只調用了 salmon quant這一步。

以下是幾個不同拼接結果的評估比較

bridger

khmer treatment and bridger

和trinity拼接結果對比

自行檢驗

可以使用 salmon 或者 kaillsto 進行表達量的快速統計分析。

至此，已經完成了常規的轉錄組拼接工作，可以進行更多的後續分析。比如基因結構註釋等等。

基因結構註釋

使用PAPS進行GENE結構註釋，一定要提前安裝好gmap或者blat中的一個或者全部

首先建立conf.txt檔案

cp $PASAHOME/pasa_conf/pasa.CONFIG.template $PASAHOME/pasa_conf/conf.txt
 ##修改conf.txt以下配置
 #MYSQLSERVER=localhost
 #MYSQL_RW_USER=mysql
 #MYSQL_RW_PASSWORD=1234
 #MYSQL_RO_USER=readonly
 #MYSQL_RO_PASSWORD=1234

配置具體任務的alignAssembly.config

cp $PASAHOME/pasa_conf/pasa.alignAssembly.Template.txt./alignAssembly.config
#修改alignAssembly.config 的內容
#MYSQLDB=<you_task_name>

用trinity生成的檔案進行基因註釋；-h檢視幫助文件。這一步使用的是適合長序列比對的軟體GMAP

/PASApipeline-2.1.0/scripts/Launch_PASA_pipeline.pl-c alignAssembly.config-C-R-g genome.fa-tTrinity.fasta--ALIGNERS gmap--transcribed_is_aligned_orient--stringent_alignment_overlap30.0--CPU5

結果檔案

{name}#.pasa_assemblies.bed/gff3/gtf/described.txt
{name}.assemblie.fasta