基於PASA進行基因預測

PASA, acronym for Program to Assemble Spliced Alignments, is a eukaryotic genome annotation tool that exploits spliced alignments of expressed transcript sequences to automatically model gene structures, and to maintain gene structure annotation consistent with the most recently available experimental sequence data. PASA also identifies and classifies all splicing variations supported by the transcript alignments.

Note:
Combine genome and Trinity de novo RNA-Seq assemblies to generate a comprehensive transcript database.

基於RNA-seq數據，將其利用trinity組裝---》利用PASA將組裝好的序列比對到draft ref----〉預測基因

1、軟件安裝

• Trinity

conda create -n trinity trinity=2.8.5
# 激活環境
conda activate trinity

•  PASA

• Mysql數據庫的准備（考慮到Mysql配置需要root權限，非常不方便，推介使用SQLite， 比起MySQL速度較慢）

1 ## 創建只讀權限（read-only）用戶和所有權限(read-write)用戶各一個
2 mysql> GRANT SELECT ON *.* TO 'pasa'@'%' IDENTIFIED BY '123456'
3 mysql> GRANT ALL ON *.* TO 'shehb'@'%' IDENTIFIED BY '123456'
4 mysql> FLUSH PRIVILEGES 123456

 

• 安裝perl模塊

cpanm install DBD::mysql
cpanm install GD
cpanm install DBD::SQLite (Sqlite需要)

 

• 安裝Gmap，blat，fasta3

conda install -c bioconda gmap
conda install blat
conda install fasta3

 

• 安裝univec 庫（污染數據庫），並建庫

wget ftp://ftp.ncbi.nih.gov/pub/UniVec/UniVec
formatdb  -i UniVec -p F

 

• 安裝PASA

## 安裝
wget https://github.com/PASApipeline/PASApipeline/releases/download/pasa-v2.4.1/PASApipeline.v2.4.1.FULL.tar.gz
tar -zxf PASApipeline.v2.4.1.FULL.tar.gz
cd PASApipeline.v2.4.1.FULL.tar
make -j 8

 

• 配置PASA config （如果使用SQLite，則下面MySQL的設置不重要）

## 配置
cd pasa_conf
cp pasa.CONFIG.template conf.txt
vi conf.txt

## 需要修改如下內容：
MYSQL_RW_USER=shehb
MYSQL_RW_PASSWORD=123456
MYSQL_RO_USER=pasa
MYSQL_RO_PASSWORD=123456
MYSQLSERVER=localhost  此處不能填寫IP
PASA_ADMIN_EMAIL=郵箱
BASE_PASA_URL=http://pasa-dev.tigr.org/cgi-bin/

 

• 修改pasa.alignAssembly.Template.txt

cd pasa_conf
cp pasa.alignAssembly.Template.txt alignAssembly.config
vi alignAssembly.config

DATABASE=/tem/mydb.sqlite
validate_alignments_in_db.dbi:--MIN_PERCENT_ALIGNED=80
validate_alignments_in_db.dbi:--MIN_AVG_PER_ID=80

　　**小提示：

• MYSQLDB指定的為工作數據庫運行完一次后需要刪除該數據庫, 否則下次運行會報錯 ；
• 如果DATABASE設置的為絕對路徑（ie. /tem/mydb.sqlite）,將會利用SQLite；如果僅僅用一個簡單名字（ie. my_pasa_db），則默認使用MySQL

 

2、運行

• 轉錄組組裝（Trinity de novo）

## 結果trinity_out_dir/Triity.fasta
Trinity --seqType fq --max_memory 50G --left reads_1.fq  --right reads_2.fq --CPU 6
# 參數
--seqType <string>      :type of reads: ('fa' or 'fq')
--max_memory <string>      :suggested max memory to use by Trinity where limiting can be enabled
--left  <string>    :left reads, one or more file names (separated by commas, no spaces)
--right <string>    :right reads, one or more file names (separated by commas, no spaces)

• 轉錄組組裝（Trinity genome-guieded)

1 Trinity --genome_guided_bam rnaseq_alignments.csorted.bam \
2         --max_memory 50G \
3         --genome_guided_max_intron 10000 \
4          --CPU 6                        
5 
6 
7 ## 說明
8 若有多個樣本的sorted.bam文件，則需使用samtools  merge 將其merge
9 結果為Trinity_GG.fasta

 

以上兩種方法可以選用一種作為transcript.fasta即可

 

• 過濾轉錄本序列（可選）

/PASApipeline.v2.4.1/bin/seqclean transcript.fasta -v /path/to/your/UniVec

 

得到transcript.fasta.cln, transcript.fasta.clean

 

• PASA將Trinity組裝結果回帖到參考基因組

Launch_PASA_pipeline.pl \
                    -c alignAssembly.config\
                    -C -R  \
                    -g example.fa.masked \
                    -t transcript.fasta.clean \
 6　　　　　　　　　　　　 -T -u transcript.fasta \
                    --ALIGNERS blat,gmap  \
                    --CPU 12 

## 參數
-c <filename> *  比對配置文件
-C  flag, create MYSQL database
-R  flag, run alignment/assembly pipeline
-g <filename> *  參考基因組
-t <filename> *  轉錄組組裝的序列
--ALIGNERS  比對軟件，可以只用一個
-- CUP 線程數

 

這一步得到的<prefix>.assemblies.fasta和<prefix>.pasa_assemblies.gff3, 其中gff3用於后面分析

 

 

====================分割線=====================

 

 此外，也可根據PASA將轉錄組回帖到基因組的結果，從中提取ORF，用於訓練基因集合，來用於其它基因預測軟件。比如：AUGUSTUS

• 從PASA組裝中提取ORF

PASApipeline-v2.3.3/scripts/pasa_asmbls_to_training_set.dbi \
    --pasa_transcripts_fasta <prefix>.assemblies.fasta \
    --pasa_transcripts_gff3 <prefix>.pasa_assemblies.gff3


##結果
<prefix>.assemblies.fasta.transdecoder.cds/pep/gff3/bed: 雖然不再基因組上，但是根據轉錄本信息，有可能是編碼區的結果
<prefix>.assemblies.fasta.transdecoder.genome.bed/gff3: 對應基因組序列的基因模型

 

 我們需要的是后者，並對其進行格式轉化gff3---> Genbank格式，進行augustus訓練，---具體可查看 Augugtus基因注釋

  

 

ref

基因結構注釋軟件PASA安裝全紀錄

PASA的安裝與使用

使用MAKER進行基因注釋(高級篇之AUGUSTUS模型訓練

Githup

 

有問題可掃描下面二維碼進行交流