EDIT und eine Lösung
Da meine ursprüngliche Frage schlecht formuliert war und ich versucht habe, das Rad neu zu erfinden, beantworte ich jetzt meine eigene Frage (vielleicht hilft es jemand anderem):
gff2fasta ist ein Tool, das genau das tut, was ich brauche, nämlich ein bestimmtes Stück DNA aus einem vollständigen Genom extrahieren (eine riesige Datei im Fasta-Format namens FULLGENOME.fasta).
Wenn ich weiß, wo sich das gewünschte Stück befindet, kann ich eine Datei mit dem Namen TEST.gff erstellen, in der ich das Gerüst (hier: sca_5_chr5_3_0) sowie den Anfang (hier: 2390621) und das Ende (hier: 2391041) meines Stücks für gff2fasta angebe, zum Beispiel:
sca_5_chr5_3_0 JGI CDS 2390621 2391041 . + 0 name "e_gw1.5.88.1"; proteinId 40463;
Dann muss ich nur noch laufen
gff2fasta.pl -gff TEST.gff -fasta FULLGENOME.fasta -out test
gff2fasta holt sich dann die Informationen aus TEST.gff und gibt den gewünschten Teil in der Datei mit dem Namen „Test“ aus, die folgendermaßen aussieht:
>sca_5_chr5_3_0_2390621_2391041_F
ATGACGACCCGCGCACCCAAAGACACATACGCTCAGCCCGACTATGAGGAGGCTCACCTAGCGACGTTTGCAGCCCCAAA
AGGCTACCCTATCGAGTCTATGCTACCCCCTAGCGTGAAGAGGGAGACCTTTGAACAGGCCCTAGCCGAGTTTACCGACG
CTATCGGCAAAGATTATGTCTTTATCGGCGATGCTCTCTCTCACTACATCGATCCCTACGACATCTATATCGATGATAGT
GAGAGAAGGAGGATGCCGAGCGCGGCTGTTTGGTACGTAACGCGGCACCCATCGAAACCTAGCAGCACTGACAAGTTTCC
GCGTCTAGTCCCTCTTCGCTCGAGGAAGCTAAGCAGGCTCTCAAGGTTGCGAACAAATACGGCATTCCGATCTGGGCATT
TTCCAGAGGCAAGAACCTGGG
Danke an Terdon und alle anderen für die Hilfe!
-------------
Dies ist die Fortsetzung dieser Frage mit weiteren Informationen und DetailsUnix: Holen Sie sich die Zeichen 10 bis 80 in einer Datei
Ich glaube, ich bin fast am Ziel, brauche aber noch etwas Hilfe.
Ich habe versucht, es so klar wie möglich zu erklären, bin mir aber bewusst, dass es sich um eine sehr spezielle Frage handelt. Lassen Sie es mich also bitte wissen, wenn ich etwas weiter klären kann!
Ich habe drei Dateien:
Hinweis an den Administrator: Ist es möglich, Dateien hochzuladen? Ich habe keine Ahnung, wie das geht ...
EinsDatei(N_haematococca_1.fasta), aus dem ich einen Namen extrahieren muss:
head -1 N_haematococca_1.fasta | cut -f4 -d'|'
dieser Name lautet in diesem Fall:
e_gw1.5.88.1
Problem 1: Der obige Code funktioniert einwandfrei, aber ich habe Probleme, den Namen (e_gw1.5.88.1) in einer Variablen zu speichern, die ich später verwenden kann ...
Ich möchte diesen Namen in einer Variablen speichern, nennen wir sie:
firstline
AzweiteDatei (Necha2_best_models.gff), in der ich alle Zeilen haben möchte, in denen dieser Name vorkommt:
grep -ir "e_gw1.5.88.1" --include="Necha2*.gff" > Necha2_in_genome.list
aber mit der benannten Variable:
grep -ir $firstline --include="Necha2*.gff" > Necha2_in_genome.list
Dies funktioniert für die Verwendung von „e_gw1.5.88.1“ …
Die Datei, die ich hier erstelle, sagt mir den Namen des DNA-Fragments, das ich schneiden möchte (sca_5_chr5_3_0) und welches Stück davon ich brauche (von Zeichennummer 2390621 bis Zeichennummer 2392655). Ich brauche diese Informationen, um die richtigen Teile aus der dritten Datei zu bekommen.
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))
In einemdritteDatei aus der ich bestimmte Teile ausschneiden möchte nach einem Schlüsselwort (sca_5_chr5_3_0) in diesem Fall:
Dank Kamaraj und hschou habe ich jetzt eine Teillösung dafür:
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
Wenn ich dies jedoch mit kleinen Zahlen debugge:
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
Ich erhalte diese Ausgabe:
r1_1_0
CCTTATCCTAGCG
nmapped
CTTATATATTAT
nmapped
TAAAAGGAGTTA
unmapped
TCTTATATAAA
unmapped
AATCTTAAGAA
Es scheint, dass die RS-Option ignoriert wird und nur die Zeichen 10 bis 20 für bestimmte Zeilen gedruckt werden. Ich habe keine Ahnung, warum diese Zeilen ausgewählt sind.
sca_5_chr5_3_0
kommt nur einmal in der Datei vor.
andere Namen, die es gibt, sind
>sca_66_unmapped
>sca_67_unmapped
usw
Ich muss diese Informationen aus 178 Genomen extrahieren. Es handelt sich dabei um riesige Dateien, und eine manuelle Suche ist einfach keine Option.
WIE die Dateien aussehen:
N_haematococca_1.fasta (Datei 1) ist eine normale Fasta-Datei:
>jgi|Necha2|40463|e_gw1.5.88.1
MTTRAPKDTYAQPDYEEAHLATFAAPKGYPIESMLPPSVKRETFEQALAEFTDAIGKDYVFIGDALSHYI
DPYDIYIDDSERRRMPSAAVCPSSLEEAKQALKVANKYGIPIWAFSRGKNLGYGGPSARVNGSVAFDLHR
Necha2_best_models.gff (Datei 2) sieht so aus (nur länger):
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 (Datei 3) sieht ungefähr so aus (viel längere GATC-Teile ...):
>sca_8_chr1_1_0
CCTTATCCTAGCGAGGATTAAGGTCCTCGAAAAGAAAGGCAAAGATATAAAGGTAATATAATAGAAGATT
AAGGTATTCTAAGTAAAGGTTATAAAGAAATAAAATAAGAAGAATATTTATAGGCTAAGAAAGACCCCCC
TAAAGGTTAAGGACTTAATATTAAGATTTAATATTCCCTAATTAATTAATATATTAATAAAAATAAAGAT
>sca_5_chr5_3_0
ATGACCACTATATCCATCGGCACAACGGCGTTAATCACATTTGGGTCTGCAATTTTCTGTTTTTGCGGTT
TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAGACTCGCTCTGCTTT
>sca_67_unmapped
CTACAGGAAGCCCTAGAGGCCCTGCAAGACCTTCCAAAGCAGCACTCTTTGCTTCTTCTTAGAGGCAGTA
TCCAGCTTCTTCTAAGGCACCTCCAGAGGCAGCTAGACCCAACCGGGCTAAAAGACCTTTGGGAAGAGGC
TAATACCCTTATAAGAGAGGCTATTATAGCCCTAGTGGCTAGAAGCCCTAGTGAGAGCCCTAAAGAGCCT
AATTCAAGCCTTATAGCCCTTCCAGTCAGAGAGGGAGGCCTAGGAATACCCTTACACAAGGACCTAGCCC
Erwartetes Endergebnis: Ist ein einzelner Text nach dem Schlüsselwort >sca_5_chr5_3_0
TTCATGTCGGTTCCAACGGGCGGAGCTCGACCAAGAATTACATACATTACGTGGCGCGAGCGCTCTACCC
TCTGGGACATGAAGGGGACGAGGAGTCTGGAGAGGTTCATGTTTGGAATCACAGCATGGTACTGCGGCAC
CCTACCTTGGTCATTCGTTGGTCTCTTTTTTTCAGGGTATGGACGCAGACTCTCAG