EDIT и решение
Поскольку мой изначальный вопрос был сформулирован неудачно и я пытался изобрести велосипед, я сейчас отвечу сам на свой вопрос (возможно, это поможет кому-то еще):
gff2fasta — это инструмент, который делает именно то, что мне нужно, а именно извлекает заданный фрагмент ДНК из полного генома (огромный файл в формате fasta, называемый FULLGENOME.fasta).
Если я знаю, где находится нужная мне часть, я могу создать файл TEST.gff, в котором укажу основу (здесь: sca_5_chr5_3_0), а также начало (здесь: 2390621) и конец (здесь: 2391041) моей части для gff2fasta, например:
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
Авторойфайл (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"...
Файл, который я создаю здесь, сообщает мне имя фрагмента ДНК, который я хочу вырезать (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), в данном случае:
Благодаря Камараджу и 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
и т. д.
Мне нужно получить эту информацию из 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