編集と解決策
私の最初の質問は表現が悪く、車輪の再発明をしようとしていたので、今は自分の質問に答えています（他の人の役に立つかもしれません）。
gff2fasta は、完全なゲノム (FULLGENOME.fasta と呼ばれる fasta 形式の巨大なファイル) から特定の DNA 部分を抽出するという、まさに私が必要としていることを実行するツールです。
必要なピースがどこにあるかわかっている場合は、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までの文字を取得する
もうすぐそこだと思うのですが、まだ助けが必要です。
できるだけわかりやすく説明しようとしましたが、非常に専門的な質問であることは承知していますので、さらに説明できることがあればお知らせください。
私が持っているのは 3 つのファイルです:
管理者へのメモ: ファイルをアップロードすることは可能ですか? 方法が全くわかりません...
1つファイル(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
あ2番この名前が出現するすべての行が必要なファイル (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 まで) が示されています。3 番目のファイルから適切な部分を取り出すには、この情報が必要です。
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
ファイル内で 1 回だけ発生します。
他にも名前はあります
>sca_66_unmapped
>sca_67_unmapped
等
この情報は 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 キーワードの後の1ビットのテキストです
TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAG