使用 gff2fasta 取代 bash 腳本從完整基因組中取得部分 DNA 序列

使用 gff2fasta 取代 bash 腳本從完整基因組中取得部分 DNA 序列

編輯和解決方案

因為我原來的問題措辭很糟糕,而且我試圖重新發明輪子,所以我現在正在回答我自己的問題(也許它對其他人有幫助):

gff2fasta 是一個完全滿足我需要的工具,即從完整基因組中提取給定的 DNA 片段(一個名為 FULLGENOME.fasta 的 fasta 格式的巨大文件)。

如果我知道我想要的片段在哪裡,我可以創建一個名為TEST.gff 的文件,在其中指定gff2fasta 片段的腳手架(此處:sca_5_chr5_3_0)以及開頭(此處:2390621)和結尾(此處: 2391041),例如:

 sca_5_chr5_3_0 JGI CDS 2390621 2391041 .   +   0   name "e_gw1.5.88.1"; proteinId 40463;

然後我只需要運行

 gff2fasta.pl -gff TEST.gff -fasta FULLGENOME.fasta -out test

然後 gff2fasta 將從 TEST.gff 獲取資訊並在名為 test 的檔案中輸出我想要的位,該檔案如下所示:

 >sca_5_chr5_3_0_2390621_2391041_F
 ATGACGACCCGCGCACCCAAAGACACATACGCTCAGCCCGACTATGAGGAGGCTCACCTAGCGACGTTTGCAGCCCCAAA
 AGGCTACCCTATCGAGTCTATGCTACCCCCTAGCGTGAAGAGGGAGACCTTTGAACAGGCCCTAGCCGAGTTTACCGACG
 CTATCGGCAAAGATTATGTCTTTATCGGCGATGCTCTCTCTCACTACATCGATCCCTACGACATCTATATCGATGATAGT
 GAGAGAAGGAGGATGCCGAGCGCGGCTGTTTGGTACGTAACGCGGCACCCATCGAAACCTAGCAGCACTGACAAGTTTCC
 GCGTCTAGTCCCTCTTCGCTCGAGGAAGCTAAGCAGGCTCTCAAGGTTGCGAACAAATACGGCATTCCGATCTGGGCATT
 TTCCAGAGGCAAGAACCTGGG

感謝 terdon 和其他人的幫助!

-------------

這是這個問題的延續,包含更多資訊和細節unix:取得文件中的第10到80個字符

我想我已經快到了,但仍然需要一些幫助。

我試圖盡可能清楚地解釋它,但我知道這是一個非常專業的問題,所以請讓我知道我是否可以進一步澄清!

我擁有的是三個文件:

管理員注意:可以上傳文件嗎?我不知道如何...

文件(N_haematococca_1.fasta) 我需要從中提取一個名稱:

 head -1 N_haematococca_1.fasta | cut -f4 -d'|' 

在本例中這個名稱是:

 e_gw1.5.88.1

問題 1:上面的程式碼工作正常,但我在將名稱 (e_gw1.5.88.1) 保存在以後可以使用的變數中時遇到問題...

我想將該名稱保存在一個變數中,我們稱之為:

 firstline

A第二檔案(Necha2_best_models.gff),我想要該名稱出現的所有行:

 grep -ir "e_gw1.5.88.1" --include="Necha2*.gff" > Necha2_in_genome.list

但對於命名變數:

 grep -ir $firstline --include="Necha2*.gff" > Necha2_in_genome.list

這適用於使用“e_gw1.5.88.1”...

我在這裡創建的文件告訴我要剪切的 DNA 片段的名稱 (sca_5_chr5_3_0) 以及我需要其中的哪一部分(從字符號 2390621 到字符號 2392655)。我需要這些資訊才能從第三個文件中獲得正確的位

 a=(`awk -F '\t' '{print $4}' Necha2_in_genome.list`)

 startDNA=${a[1]}
 endDNA=${a[${#a[@]}-1]}

 #add 1000 or other number, depending on the problems with the gene
 correctedstartDNA=$(($startDNA-1000))
 correctedendDNA=$(($endDNA-1000))

在一個第三在本例中,我想從中剪下關鍵字(sca_5_chr5_3_0)後的某些部分:

感謝 Kamaraj 和 hschou,我現在對此有了部分解決方案:

 cat Necha2_scaffolds.fasta | sed -n -e '/sca_5_chr5_3_0/,$p' | grep -v '^>' | tr -d '\n'|awk -v start="${correctedstartDNA}" -v end="${correctedendDNA}" '{print substr($0,start,end)}' RS= Necha2_scaffolds.fasta

但是,如果我用較小的數字進行調試:

 cat Necha2_scaffolds.fasta | sed -n -e '/sca_5_chr5_3_0/,$p' | grep -v '^>' |  tr -d '\n'|awk -v start=10 -v end=20 '{print substr($0,start,end)}' RS= Necha2_scaffolds.fasta

我得到這個輸出:

 r1_1_0
 CCTTATCCTAGCG
 nmapped
 CTTATATATTAT
 nmapped
 TAAAAGGAGTTA
 unmapped
 TCTTATATAAA
 unmapped
 AATCTTAAGAA

看來 RS 選項被忽略,它只列印某些行的字元 10 到 20。我不知道為什麼選擇這些行。

 sca_5_chr5_3_0

僅在文件中出現一次。

還有其他名字

 >sca_66_unmapped 
 >sca_67_unmapped

ETC

我必須從 178 個基因組中獲取這些信息,它們都是巨大的文件,手動搜索不是一個選擇。

文件看起來如何:

N_haematococca_1.fasta(文件1)是一個普通的fasta檔案:

 >jgi|Necha2|40463|e_gw1.5.88.1
 MTTRAPKDTYAQPDYEEAHLATFAAPKGYPIESMLPPSVKRETFEQALAEFTDAIGKDYVFIGDALSHYI
 DPYDIYIDDSERRRMPSAAVCPSSLEEAKQALKVANKYGIPIWAFSRGKNLGYGGPSARVNGSVAFDLHR

Necha2_best_models.gff(檔案 2)看起來像這樣(只是更長):

 sca_5_chr5_3_0 JGI exon    2390621 2390892 .   +   .   name "e_gw1.5.88.1"; transcriptId 40463
 sca_5_chr5_3_0 JGI CDS 2390621 2390892 .   +   0   name "e_gw1.5.88.1"; proteinId 40463; exonNumber 1
 sca_5_chr5_3_0 JGI start_codon 2390621 2390623 .   +   0   name "e_gw1.5.88.1"
 sca_5_chr5_3_0 JGI exon    2390949 2391041 .   +   .   name "e_gw1.5.88.1"; transcriptId 40463
 sca_5_chr5_3_0 JGI CDS 2390949 2391041 .   +   2   name "e_gw1.5.88.1"; proteinId 40463; exonNumber 2
 sca_5_chr5_3_0 JGI exon    2391104 2391311 .   +   .   name "e_gw1.5.88.1"; transcriptId 40463
 sca_5_chr5_3_0 JGI CDS 2391104 2391311 .   +   2   name "e_gw1.5.88.1"; proteinId 40463; exonNumber 3
 sca_5_chr5_3_0 JGI exon    2391380 2392367 .   +   .   name "e_gw1.5.88.1"; transcriptId 40463
 sca_5_chr5_3_0 JGI CDS 2391380 2392367 .   +   0   name "e_gw1.5.88.1"; proteinId 40463; exonNumber 4
 sca_5_chr5_3_0 JGI exon    2392421 2392485 .   +   .   name "e_gw1.5.88.1"; transcriptId 40463
 sca_5_chr5_3_0 JGI CDS 2392421 2392485 .   +   1   name "e_gw1.5.88.1"; proteinId 40463; exonNumber 5
 sca_5_chr5_3_0 JGI exon    2392541 2392657 .   +   .   name "e_gw1.5.88.1"; transcriptId 40463
 sca_5_chr5_3_0 JGI CDS 2392541 2392657 .   +   0   name "e_gw1.5.88.1"; proteinId 40463; exonNumber 6
 sca_5_chr5_3_0 JGI stop_codon  2392655 2392657 .   +   0   name "e_gw1.5.88.1"
 sca_5_chr5_3_0 JGI exon    2396205 2396730 .   +   .   name "e_gw1.5.227.1"; transcriptId 41333
 sca_5_chr5_3_0 JGI CDS 2396205 2396730 .   +   0   name "e_gw1.5.227.1"; proteinId 41333; exonNumber 1
 sca_5_chr5_3_0 JGI start_codon 2396205 2396207 .   +   0   name "e_gw1.5.227.1"

Necha2_scaffolds.fasta(文件 3)看起來有點像這樣(更長的 GATC 位元...):

 >sca_8_chr1_1_0
 CCTTATCCTAGCGAGGATTAAGGTCCTCGAAAAGAAAGGCAAAGATATAAAGGTAATATAATAGAAGATT
 AAGGTATTCTAAGTAAAGGTTATAAAGAAATAAAATAAGAAGAATATTTATAGGCTAAGAAAGACCCCCC
 TAAAGGTTAAGGACTTAATATTAAGATTTAATATTCCCTAATTAATTAATATATTAATAAAAATAAAGAT
 >sca_5_chr5_3_0
 ATGACCACTATATCCATCGGCACAACGGCGTTAATCACATTTGGGTCTGCAATTTTCTGTTTTTGCGGTT
 TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
 TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
 CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAGACTCGCTCTGCTTT
 >sca_67_unmapped
 CTACAGGAAGCCCTAGAGGCCCTGCAAGACCTTCCAAAGCAGCACTCTTTGCTTCTTCTTAGAGGCAGTA
 TCCAGCTTCTTCTAAGGCACCTCCAGAGGCAGCTAGACCCAACCGGGCTAAAAGACCTTTGGGAAGAGGC
 TAATACCCTTATAAGAGAGGCTATTATAGCCCTAGTGGCTAGAAGCCCTAGTGAGAGCCCTAAAGAGCCT
 AATTCAAGCCTTATAGCCCTTCCAGTCAGAGAGGGAGGCCTAGGAATACCCTTACACAAGGACCTAGCCC

預期最終輸出:是 >sca_5_chr5_3_0 關鍵字之後的一位文本

 TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
 TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
 CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAG

相關內容