轉(zhuǎn)錄本預(yù)測編碼區(qū)：在transcript中找CDS區(qū)

RNA-Seq數(shù)據(jù)組裝得到的轉(zhuǎn)錄本具有編碼蛋白質(zhì)的能力咱圆。預(yù)測編碼區(qū)對于確定轉(zhuǎn)錄產(chǎn)物在細(xì)胞中的分子作用至關(guān)重要!!!

通過Cufflinks或Stringtie得到的gtf文件也物，所以需要準(zhǔn)備如下兩個(gè)文件：

  1. transcripts.gtf: 你自己的RNA-seq數(shù)據(jù)經(jīng)過Cufflinks控淡，或Stringtie得到的gtf文件陪白，記錄預(yù)測轉(zhuǎn)錄本的GTF文件

  2. genome.fasta: 參考基因組序列

第一步：依據(jù)參考genome.fasta文件，從gtf文件中提取fasta序列叠荠。

conda activate python3 #進(jìn)入環(huán)境
TransDecoder.Predict #啟動TransDecoder

###基于mm10的 Tophat結(jié)果
示例：
gtf_genome_to_cdna_fasta.pl trans.gtf /xx/xx/Mus_musculus.GRCm38.dna.toplevel.fa > trans.fa
操作路徑于：/home/u20111230014/workspace/tran_gtf_ZF 
PS: DPP-0_transcripts.gtf和 mm10.fa必須在同一個(gè)文件夾
實(shí)例：
gtf_genome_to_cdna_fasta.pl DPP-0_transcripts.gtf mm10.fa >  DPP-0_mm10_trans.fa
gtf_genome_to_cdna_fasta.pl DPP-1_transcripts.gtf mm10.fa >  DPP-1_mm10_trans.fa

###基于GRCm39的 Tophat結(jié)果
操作所在路徑：/home/u20111230014/workspace/genome/GRCm39/tophat_gtf_genome_to_cdna

實(shí)例
將所有所需要的文件軟連接到同一文件夾
ln -s /home/u20111230014/workspace/genome/GRCm39/Tophat_Cufflinks_DPP-0_clean.fq/transcripts.gtf DPP-0_GRCm39_transcripts.gtf 

gtf_genome_to_cdna_fasta.pl DPP-0_GRCm39_transcripts.gtf GRCm39.fa >  DPP-0_GRCm39_trans.fa
gtf_genome_to_cdna_fasta.pl DPP-1_GRCm39_transcripts.gtf GRCm39.fa >  DPP-1_GRCm39_trans.fa

第二步：將gtf文件轉(zhuǎn)換為gff3格式(只是為了后續(xù)分析)盹沈。

示例
gtf_to_alignment_gff3.pl trans.gtf > trans.gff3

實(shí)例
###基于mm10的 Tophat結(jié)果
操作路徑于：/home/u20111230014/workspace/tran_gtf_ZF 
gtf_to_alignment_gff3.pl DPP-0_transcripts.gtf >  DPP-0_mm10_trans.gff3
gtf_to_alignment_gff3.pl DPP-1_transcripts.gtf >  DPP-1_mm10_trans.gff3

###基于GRCm39的 Tophat結(jié)果
操作所在路徑：/home/u20111230014/workspace/genome/GRCm39/tophat_gtf_genome_to_cdna
軟連接腳本：

ln -s /home/u20111230014/workspace/tran_gtf_ZF/gtf_to_alignment_gff3.slurm tophat_gtf_to_alignment_gff3.slurm
實(shí)例
gtf_to_alignment_gff3.pl DPP-0_GRCm39_transcripts.gtf >  DPP-0_GRCm39_trans.gff3
gtf_to_alignment_gff3.pl DPP-1_GRCm39_transcripts.gtf >  DPP-1_GRCm39_trans.gff3

第三步：預(yù)測轉(zhuǎn)錄本中長的開放閱讀框，-m參數(shù)設(shè)置為300個(gè)氨基酸袜香。并進(jìn)行blastp和hmmer比對撕予。

示例
TransDecoder.LongOrfs -m 300 -t trans.fa
實(shí)例
TransDecoder.LongOrfs -m 300 -t DPP-0_mm10_trans.fa
TransDecoder.LongOrfs -m 300 -t DPP-1_mm10_trans.fa

第四步：為了提高預(yù)測的靈敏度，與已知的蛋白數(shù)據(jù)庫進(jìn)行比對蜈首，探索同源性实抡。

1) blastp search

選用Swissprot數(shù)據(jù)庫的蛋白序列，下載uniprot_sprot.fasta.gz欢策，網(wǎng)址是ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/

解壓縮

 gunzip uniprot_sprot.fasta.gz

建blast+庫于路徑：/home/u20111230014/workspace/genome  
#### 建數(shù)據(jù)庫
makeblastdb -in uniprot_sprot.fasta -dbtype prot -parse_seqids -hash_index -out uniprot

## makeblastdb參數(shù)說明：

-in 輸入數(shù)據(jù)庫文件

-dbtype 蛋白庫用prot,核酸用nucl

-out 所建數(shù)據(jù)庫的名稱

 -parse_seqids      

-hash_index 創(chuàng)建哈希序列

數(shù)據(jù)比對

###基于mm10.fa的比對
操作路徑在：/home/u20111230014/workspace/tran_gtf_ZF 
index文件：uniprot的11個(gè)子文件吆寨，再加上一個(gè)uniprot_sprot.fasta，共12個(gè)
示例：
blastp -query trans.fa.transdecoder_dir/longest_orfs.pep -db uniprot -max_target_seqs 1 -outfmt 6 -evalue 1e-5 -num_threads 10 > blastp.outfmt6 &
實(shí)例：
blastp -query DPP-0_mm10_trans.fa.transdecoder_dir/longest_orfs.pep -db uniprot -max_target_seqs 1 -outfmt 6 -evalue 1e-5 -num_threads 10 > DPP-0_mm10_blastp.outfmt6

２瓤堋Ｗ那濉！PS: 如果報(bào)錯(cuò)：BLAST Database error: No alias or index file found for nucleotide database [nt] in search path, 必須將uniprot的index俺孙，共12個(gè)文件辣卒，跟要跑的所有文件放在一起才行，如下：

(python3) [u20111230014@cpu10 tran_gtf_ZF]$ ll uniprot*
lrwxrwxrwx 1 u20111230014 u20111230014 43 Feb 23 16:23 uniprot -> /home/u20111230014/workspace/genome/uniprot
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.pdb -> /home/u20111230014/workspace/genome/uniprot/uniprot.pdb
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.phd -> /home/u20111230014/workspace/genome/uniprot/uniprot.phd
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.phi -> /home/u20111230014/workspace/genome/uniprot/uniprot.phi
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.phr -> /home/u20111230014/workspace/genome/uniprot/uniprot.phr
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.pin -> /home/u20111230014/workspace/genome/uniprot/uniprot.pin
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.pog -> /home/u20111230014/workspace/genome/uniprot/uniprot.pog
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.pos -> /home/u20111230014/workspace/genome/uniprot/uniprot.pos
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.pot -> /home/u20111230014/workspace/genome/uniprot/uniprot.pot
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.psq -> /home/u20111230014/workspace/genome/uniprot/uniprot.psq
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.ptf -> /home/u20111230014/workspace/genome/uniprot/uniprot.ptf
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.pto -> /home/u20111230014/workspace/genome/uniprot/uniprot.pto
lrwxrwxrwx 1 u20111230014 u20111230014 63 Feb 23 16:28 uniprot_sprot.fasta -> /home/u20111230014/workspace/genome/uniprot/uniprot_sprot.fasta

如果還是報(bào)同樣的錯(cuò)睛榄，就換一種方法：下載preformatted的數(shù)據(jù)庫添寺！

去到miniconda下 blast包里含update_blastdb.pl文件的路徑，用perl 更新一下數(shù)據(jù)庫
!!!blast軟件所在的位置：/home/u20111230014/miniconda3/pkgs/blast-2.12.0-h3289130_3/bin

先查看數(shù)據(jù)庫列表：
(python3) [u20111230014@cpu10 bin]$ perl update_blastdb.pl --showall
Connected to NCBI
16S_ribosomal_RNA
18S_fungal_sequences
28S_fungal_sequences
Betacoronavirus
ITS_RefSeq_Fungi
ITS_eukaryote_sequences
LSU_eukaryote_rRNA
LSU_prokaryote_rRNA
SSU_eukaryote_rRNA
env_nt
env_nr
human_genome
landmark
mito
mouse_genome
nr
nt
pataa
patnt
pdbaa
pdbnt
ref_euk_rep_genomes
ref_prok_rep_genomes
ref_viroids_rep_genomes
ref_viruses_rep_genomes
refseq_select_rna
refseq_select_prot
refseq_protein
refseq_rna
swissprot
tsa_nr
tsa_nt
taxdb
###再選擇需要的數(shù)據(jù)庫進(jìn)行更新下載：
perl update_blastdb.pl nr #nr 是蛋白質(zhì)懈费， nt是核酸
!!!PS:nohup掛在后臺下載计露，中途如果超時(shí)的話，需重新執(zhí)行命令就行，之前下載好的不用重新執(zhí)行Ｆ惫蕖２嫒ぁ！大概有十幾個(gè)文件该押，下好批量解壓命令：ls *.tar.gz | xargs -n1 tar xzvf

基于GRCm39.fa的比對

操作路徑在：/home/u20111230014/workspace/genome/GRCm39/tophat_gtf_genome_to_cdna/
index文件：uniprot的11個(gè)子文件疗杉，再加上一個(gè)uniprot_sprot.fasta，共12個(gè)
實(shí)例
blastp -query DPP-0_mm10_trans.fa.transdecoder_dir/longest_orfs.pep -db uniprot -max_target_seqs 1 -outfmt 6 -evalue 1e-5 -num_threads 10 > DPP-0_mm10_blastp.outfmt6

2) Pfam search

在建立索引數(shù)據(jù)庫hmmpress Pfam-A.hmm成功后蚕礼，使用 hmmscan 進(jìn)行 Pfam 注釋
conda activate python3 #進(jìn)入環(huán)境
hmmbuild --help #啟動
Pfam-A存放于：hmm=/home/u20111230014/workspace/genome/Pfam-A/Pfam-A.hmm
PS: 執(zhí)行hmmscan 前輸入上面的路徑指示烟具，或直接寫入slurm腳本

進(jìn)入存放轉(zhuǎn)換格式好的樣本tran.fa文件目錄

cd /home/u20111230014/workspace/tran_gtf_ZF 
(python3) [u20111230014@cpu6 tran_gtf_ZF]$ ll *mm10_trans.fa
-rw-rw-r-- 1 u20111230014 u20111230014 38711086 Feb 23 15:26 DPP-0_mm10_trans.fa
-rw-rw-r-- 1 u20111230014 u20111230014 38539068 Feb 23 15:28 DPP-1_mm10_trans.fa
-rw-rw-r-- 1 u20111230014 u20111230014 45990886 Feb 23 15:31 DPP-2_mm10_trans.fa
-rw-rw-r-- 1 u20111230014 u20111230014 39814002 Feb 23 15:34 DPP-3_mm10_trans.fa
-rw-rw-r-- 1 u20111230014 u20111230014 33385081 Feb 23 15:37 DPP-4_mm10_trans.fa
-rw-rw-r-- 1 u20111230014 u20111230014 38240368 Feb 23 15:40 DPP-5_mm10_trans.fa

示例：
hmmscan --cpu 8 -o hmm_pfam.trans.txt --tblout hmm_pfam.trans.tbl --noali -E 1e-5 $hmm trans.fa

###參數(shù)說明
 -o FILE 將結(jié)果輸出到指定的文件中。默認(rèn)是輸出到標(biāo)準(zhǔn)輸出

 --tblout FILE  將蛋白質(zhì)家族的結(jié)果以表格形式輸出到指定的文件中奠蹬。默認(rèn)不輸出該文件

--domtblout FILE  將蛋白結(jié)構(gòu)域的比對結(jié)果以表格形式輸出到指定的文件中朝聋。默認(rèn)不輸出該文件。該表格中包含query序列起始結(jié)束位點(diǎn)與目標(biāo)序列起始結(jié)束位點(diǎn)的匹配信息

 --acc 在輸出結(jié)果中包含 PF 的編號囤躁，默認(rèn)是蛋白質(zhì)家族的名稱

--noali 在輸出結(jié)果中不包含比對信息冀痕。輸出文件的大小則會更小

-E FLOAT    default:10.0 設(shè)定 E_value 閾值，推薦設(shè)置為 1e-5

 -T FLOAT  設(shè)定 Score 閾值

--domE FLOAT    default:10.0  設(shè)定 E_value 閾值狸演。該參數(shù)和 -E 參數(shù)類似言蛇，不過是 domain 比對設(shè)定的值

--cpu 線程
實(shí)例腳本：
hmm=/home/u20111230014/workspace/genome/Pfam-A/Pfam-A.hmm
hmmscan --cpu 25 -o DPP-0_mm10_hmm_pfam.trans.txt --tblout DPP-0_mm10_hmm_pfam.trans.tbl --noali -E 1e-5 $hmm DPP-0_mm10_trans.fa
hmmscan --cpu 25 -o DPP-1_mm10_hmm_pfam.trans.txt --tblout DPP-1_mm10_hmm_pfam.trans.tbl --noali -E 1e-5 $hmm DPP-1_mm10_trans.fa
####
hmmscan --cpu 25 -o DPP-0_GRCm39_hmm_pfam.trans.txt --tblout DPP-0_GRCm39_hmm_pfam.trans.tbl --noali -E 1e-5 $hmm DPP-0_GRCm39_trans.fa

hmmscan --cpu 25 -o DPP-1_GRCm39_hmm_pfam.trans.txt --tblout DPP-1_GRCm39_hmm_pfam.trans.tbl --noali -E 1e-5 $hmm DPP-1_GRCm39_trans.fa

第五步：結(jié)合blastp和pfam的結(jié)果預(yù)測可能的編碼區(qū)。

TransDecoder.Predict #啟動
示例：
TransDecoder.Predict -t trans.fa --retain_pfam_hits hmm_pfam.trans.tbl --retain_blastp_hits blastp.outfmt6
實(shí)例：
TransDecoder.Predict -t DPP-0_GRCm39_trans.fa --retain_pfam_hits DPP-0_GRCm39_hmm_pfam.trans.tbl --retain_blastp_hits DPP-0_GRCm39_blastp.outfmt6

第六步: 生成基于參考基因組的編碼區(qū)注釋文件

示例：
cdna_alignment_orf_to_genome_orf.pl trans.fa.transdecoder.gff3 trans.gff3 trans.fa > trans.fa.transdecoder.genome.gff3

實(shí)例：
cdna_alignment_orf_to_genome_orf.pl DPP-0_mm10_trans.fa.transdecoder.gff3 DPP-0_mm10_trans.gff3 DPP-0_mm10_trans.fa > DPP-0_mm10_trans.fa.transdecoder.genome.gff3

最終輸出結(jié)果文件如下:

  1. trans.fa.transdecoder.pep: 最終預(yù)測的CDS對應(yīng)的蛋白序列

  2. trans.fa.transdecoder.cds: 最終預(yù)測的CDS序列

  3. trans.fa.transdecoder.gff3: 最終ORF對應(yīng)的GFF3

  4. trans.fa.transdecoder.bed: 以BED格式存放ORF位置信息

  5. trans.fa.transdecoder.genome.gff3: 基于參考基因組的GFF3文件

BED和GFF3文件可以使用IGV可視化宵距。

參考：（[小張聊科研]）https://mp.weixin.qq.com/s?src=11&timestamp=1645587577&ver=3637&signature=ydPEWFG06TSywWPWpAOTRgJNtSQ2BUEOlUYjQZV3OGiZKZxJmIp2mHe9bWlV3UDcnKn0p1K0SWstUxzPQyTmCKC7GdYDBRdHKhKYct*vh-mH-qsE4a4A-TBNazUb&new=1

生信分析師高端
http://dl2.kerwin.cn:8063/csdn/key/article-a_giant_pig-103065967/auth/1645608443-2022062317272318-0-2af5159d4f371312acc26ac982c48187

卡西莫多的禮物
https://blog.csdn.net/qq_35696312/article/details/104779347

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末腊尚，一起剝皮案震驚了整個(gè)濱河市，隨后出現(xiàn)的幾起案子满哪，更是在濱河造成了極大的恐慌跟伏，老刑警劉巖，帶你破解...
    沈念sama閱讀 206,839評論 6贊 482
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件翩瓜，死亡現(xiàn)場離奇詭異受扳，居然都是意外死亡，警方通過查閱死者的電腦和手機(jī)兔跌，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,543評論 2贊 382
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門勘高，熙熙樓的掌柜王于貴愁眉苦臉地迎上來，“玉大人坟桅，你說我怎么就攤上這事华望。” “怎么了仅乓？”我有些...
    開封第一講書人閱讀 153,116評論 0贊 344
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵赖舟，是天一觀的道長。 經(jīng)常有香客問我夸楣，道長宾抓，這世上最難降的妖魔是什么子漩？ 我笑而不...
    開封第一講書人閱讀 55,371評論 1贊 279
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任，我火速辦了婚禮石洗，結(jié)果婚禮上幢泼，老公的妹妹穿的比我還像新娘。我一直安慰自己讲衫，他們只是感情好缕棵，可當(dāng)我...
    茶點(diǎn)故事閱讀 64,384評論 5贊 374
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布。 她就那樣靜靜地躺著涉兽，像睡著了一般招驴。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上枷畏，一...
    開封第一講書人閱讀 49,111評論 1贊 285
  • 城市分裂傳說
    那天别厘，我揣著相機(jī)與錄音，去河邊找鬼矿辽。 笑死，一個(gè)胖子當(dāng)著我的面吹牛郭厌，可吹牛的內(nèi)容都是我干的袋倔。 我是一名探鬼主播，決...
    沈念sama閱讀 38,416評論 3贊 400
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼折柠，長吁一口氣：“原來是場噩夢啊……” “哼宾娜！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起扇售，我...
    開封第一講書人閱讀 37,053評論 0贊 259
  • 萬榮殺人案實(shí)錄
    序言：老撾萬榮一對情侶失蹤前塔，失蹤者是張志新（化名）和其女友劉穎，沒想到半個(gè)月后承冰，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體华弓，經(jīng)...
    沈念sama閱讀 43,558評論 1贊 300
  • ?護(hù)林員之死
    正文 獨(dú)居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 36,007評論 2贊 325
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年困乒，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了寂屏。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點(diǎn)故事閱讀 38,117評論 1贊 334
  • 活死人
    序言：一個(gè)原本活蹦亂跳的男人離奇死亡娜搂，死狀恐怖迁霎，靈堂內(nèi)的尸體忽然破棺而出，到底是詐尸還是另有隱情百宇，我是刑警寧澤考廉，帶...
    沈念sama閱讀 33,756評論 4贊 324
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布，位于F島的核電站携御，受9級特大地震影響昌粤，放射性物質(zhì)發(fā)生泄漏既绕。R本人自食惡果不足惜，卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 39,324評論 3贊 307
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一婚苹、第九天 我趴在偏房一處隱蔽的房頂上張望岸更。 院中可真熱鬧，春花似錦膊升、人聲如沸怎炊。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,315評論 0贊 19
  • 一樁弒父案廓译，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽评肆。三九已至，卻和暖如春非区，著一層夾襖步出監(jiān)牢的瞬間瓜挽，已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 31,539評論 1贊 262
  • 情欲美人皮
    我被黑心中介騙來泰國打工征绸， 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留久橙，地道東北人。 一個(gè)月前我還...
    沈念sama閱讀 45,578評論 2贊 355
  • 代替公主和親
    正文 我出身青樓管怠，卻偏偏與公主長得像淆衷，于是被迫代替她去往敵國和親。 傳聞我的和親對象是個(gè)殘疾皇子渤弛，可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 42,877評論 2贊 345

推薦閱讀更多精彩內(nèi)容