![使用 gff2fasta 取代 bash 腳本從完整基因組中取得部分 DNA 序列](https://rvso.com/image/109212/%E4%BD%BF%E7%94%A8%20gff2fasta%20%E5%8F%96%E4%BB%A3%20bash%20%E8%85%B3%E6%9C%AC%E5%BE%9E%E5%AE%8C%E6%95%B4%E5%9F%BA%E5%9B%A0%E7%B5%84%E4%B8%AD%E5%8F%96%E5%BE%97%E9%83%A8%E5%88%86%20DNA%20%E5%BA%8F%E5%88%97.png)
編輯和解決方案
因為我原來的問題措辭很糟糕,而且我試圖重新發明輪子,所以我現在正在回答我自己的問題(也許它對其他人有幫助):
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