Скрипт с использованием sed и grep дает непреднамеренный вывод

Question

У меня есть файл «source.fasta», который содержит информацию в следующем формате:

>TRINITY_DN80_c0_g1_i1 len=723 path=[700:0-350 1417:351-368 1045:369-722] [-1, 700, 1417, 1045, -2]
CGTGGATAACACATAAGTCACTGTAATTTAAAAACTGTAGGACTTAGATCTCCTTTCTATATTTTTCTGATAACATATGGAACCCTGCCGATCATCCGATTTGTAATATACTTAACTGCTGGATAACTAGCCAAAAGTCATCAGGTTATTATATTCAATAAAATGTAACTTGCCGTAAGTAACAGAGGTCATATGTTCCTGTTCGTCACTCTGTAGTTACAAATTATGACACGTGTGCGCTG
>TRINITY_DN83_c0_g1_i1 len=371 path=[1:0-173 152:174-370] [-1, 1, 152, -2]
GTTGTAAACTTGTATACAATTGGGTTCATTAAAGTGTGCACATTATTTCATAGTTGATTTGATTATTCCGAGTGACCTATTTCGTCACTCGATGTTTAAAGAAATTGCTAGTGTGACCCCAATTGCGTCAGACCAAAGATTGAATCTAGACATTAATTTCCTTTTGTATTTGTATCGAGTAAGTTTACAGTCGTAAATAAAGAATCTGCCTTGAACAAACCTTATTCCTTTTTATTCTAAAAGAGGCCTTTGCGTAGTAGTAACATAGTACAAATTGGTTTATTTAACGATTTATAAACGATATCCTTCCTACAGTCGGGTGAAAAGAAAGTATTCGAAATTAGAATGGTTCCTCATATTACACGTTGCTG
>TRINITY_DN83_c0_g1_i2 len=218 path=[1:0-173 741:174-217] [-1, 1, 741, -2]
GTTGTAAACTTGTATACAATTGGGTTCATTAAAGTGTGCACATTATTTCATAGTTGATTTGATTATTCCGAGTGACCTATTTCGTCACTCGATGTTTAAAGAAATTGCTAGTGTGACCCCAATTGCGTCAGACCAAAGATTGAATCTAGACATTAATTTCCTTTTGTATTTGTACCGAGTAAGTTTCCAGTCGTAAATAAAGAATCTGCCAGATCGGA
>TRINITY_DN99_c0_g1_i1 len=326 path=[1:0-242 221:243-243 652:244-267 246:268-325] [-1, 1, 221, 652, 246, -2]
ATCGGTACTATCATGTCATATATCTAGAAATAATACCTACGAATGTTATAAGAATTTCATAACATGATATAACGATCATACATCATGGCCTTTCGAAGAAAATGGCGCATTTACGTTTAATAATTCCGCGAAAGTCAAGGCAAATACAGACCTAATGCGAAATTGAAAAGAAAATCCGAATATCAGAAACAGAACCCAGAACCAATATGCTCAGCAGTTGCTTTGTAGCCAATAAACTCAACTAGAAATTGCTTATCTTTTATGTAACGCCATAAAACGTAATACCGATAACAGACTAAGCACACATATGTAAATTACCTGCTAAA
>TRINITY_DN90_c0_g1_i1 len=1240 path=[1970:0-527 753:528-1239] [-1, 1970, 753, -2]
GTCGATACTAGACAACGAATAATTGTGTCTATTTTTAAAAATAATTCCTTTTGTAAGCAGATTTTTTTTTTCATGCATGTTTCGAGTAAATTGGATTACGCATTCCACGTAACATCGTAAATGTAACCACATTGTTGTAACATACGGTATTTTTTCTGACAACGGACTCGATTGTAAGCAACTTTGTAACATTATAATCCTATGAGTATGACATTCTTAATAATAGCAACAGGGATAAAAATAAAACTACATTGTTTCATTCAACTCGTAAGTGTTTATTTAAAATTATTATTAAACACTATTGTAATAAAGTTTATATTCCTTTGTCAGTGGTAGACACATAAACAGTTTTCGAGTTCACTGTCG
>TRINITY_DN84_c0_g1_i1 len=301 path=[1:0-220 358:221-300] [-1, 1, 358, -2]
ACTATTATGTAGTACCTACATTAGAAACAACTGACCCAAGACAGGAGAAGTCATTGGATGATTTTCCCCATTAAAAAAAGACAACCTTTTAAGTAAGCATACTCCAAATTAAGGTTTAATTAGCTAAGTGAGCGCGAAAAATGATCAAATATACCGACGTCCATTTGGGGCCTATCCTTTTTAGTGTTCCTAATTGAAATCCTCACGTATACAGCTAGTCACTTTTAAATCTTATAAACATGTGATCCGTCTGCTCATTTGGACGTTACTGCCCAAAGTTGGTACATGTTTCGTACTCACG
>TRINITY_DN84_c0_g1_i2 len=301 path=[1:0-220 199:221-300] [-1, 1, 199, -2]
ACTATTATGTAGTACCTACATTAGAAACAACTGACCCAAGACAGGAGAAGTCATTGGATGATTTTCCCCATTAAAAAAAGACAACCTTTTAAGTAAGCATACTCCAAATTAAGGTTTAATTAGCTAAGTGAGCGCGAAAAATGATCAAATATACCGACGTCCATTTGGGGCCTATCCTTTTTAGTGTTCCTAATTGAAATCCTCACGTATACAGCTAGTCAGCTAACCAAAGATAAGTGTCTTGGCTTGGTATCTACAGATCTCTTTTCGTAATTTCGTGAGTACGAAACATGTACCAACT
>TRINITY_DN72_c0_g1_i1 len=434 path=[412:0-247 847:248-271 661:272-433] [-1, 412, 847, 661, -2]
GTTAATTTAGTGGGAAGTATGTGTTAAAATTAGTAAATTAGGTGTTGGTGTGTTTTTAATATGAATCCGGAAGTGTTTTGTTAGGTTACAAGGGTACGGAATTGTAATAATAGAAATCGGTATCCTTGAGACCAATGTTATCGCATTCGATGCAAGAATAGATTGGGAAATAGTCCGGTTATCAATTACTTAAAGATTTCTATCTTGAAAACTATTTCTAATTGGTAAAAAAACTTATTTAGAATCACCCATAGTTGGAAGTTTAAGATTTGAGACATCTTAAATTTTTGGTAGGTAATTTTAAGATTCTATCGTAGTTAGTACCTTTCGTTCTTCTTATTTTATTTGTAAAATATATTACATTTAGTACGAGTATTGTATTTCCAATATTCAGTCTAATTAGAATTGCAAAATTACTGAACACTCAATCATAA
>TRINITY_DN75_c0_g1_i1 len=478 path=[456:0-477] [-1, 456, -2]
CGAGCACATCAGGCCAGGGTTCCCCAAGTGCTCGAGTTTCGTAACCAAACAACCATCTTCTGGTCCGACCACCAGTCACATGATCAGCTGTGGCGCTCAGTATACGAGCACAGATTGCAACAGCCACCAAATGAGAGAGGAAAGTCATCCACATTGCCATGAAATCTGCGAAAGAGCGTAAATTGCGAGTAGCATGACCGCAGGTACGGCGCAGTAGCTGGAGTTGGCAGCGGCTAGGGGTGCCAGGAGGAGTGCTCCAAGGGTCCATCGTGCTCCACATGCCTCCCCGCCGCTGAACGCGCTCAGAGCCTTGCTCATCTTGCTACGCTCGCTCCGTTCAGTCATCTTCGTGTCTCATCGTCGCAGCGCGTAGTATTTACG

В этом файле содержится около 400 000 последовательностей.

У меня есть другие идентификаторы файлов.txt в следующем формате:

>TRINITY_DN14840_c10_g1_i1
>TRINITY_DN8506_c0_g1_i1
>TRINITY_DN12276_c0_g2_i1
>TRINITY_DN15434_c5_g3_i1
>TRINITY_DN9323_c8_g3_i5
>TRINITY_DN11957_c1_g7_i1
>TRINITY_DN15373_c1_g1_i1
>TRINITY_DN22913_c0_g1_i1
>TRINITY_DN13029_c4_g5_i1

У меня есть 100 идентификаторов последовательности в этом файле.Когда я сопоставляю эти идентификаторы с исходным файлом, я хочу вывод, который дает мне совпадение для каждого из этих идентификаторов со всей последовательностью.

Например, для идентификатора:

>TRINITY_DN80_c0_g1_i1

Я хочу, чтобы мой вывод был:

>TRINITY_DN80_c0_g1_i1
CGTGGATAACACATAAGTCACTGTAATTTAAAAACTGTAGGACTTAGATCTCCTTTCTATATTTTTCTGATAACATATGGAACCCTGCCGATCATCCGATTTGTAATATACTTAACTGCTGGATAACTAGCCAAAAGTCATCAGGTTATTATATTCAATAAAATGTAACTTGCCGTAAGTAACAGAGGTCATATGTTCCTGTTCGTCACTCTGTAGTTACAAATTATGACACGTGTGCGCTG

Я хочу, чтобы все сто последовательностей в этом формате.Я использовал этот код:

while read p; do
echo ''$p >> out.fasta
grep -A 400000 -w $p source.fasta | sed -n -e '1,/>/ {/>/ !{'p''}} >> out.fasta
done < ids.txt

Но мой вывод отличается тем, что только последний идентификатор имеет последовательность, а остальные не имеют никакой связанной последовательности:

>TRINITY_DN14840_c10_g1_i1
>TRINITY_DN8506_c0_g1_i1
>TRINITY_DN12276_c0_g2_i1
....
>TRINITY_DN10309_c6_g3_i1
>TRINITY_DN6990_c0_g1_i1
TTTTTTTTTTTTTGTGGAAAAACATTGATTTTATTGAATTGTAAACTTAAAATTAGATTGGCTGCACATCTTAGATTTTGTTGAAAGCAGCAATATCAACAGACTGGACGAAGTCTTCGAATTCCTGGATTTTTTCAGTCAAGAGATCAACAGACACTTTGTCGTCTTCAATGACACACATGATCTGCAGTTTGTTGATACCATATCCAACAGGTACAAGTTTGGAAGCTCCCCAGAGGAGACCATCCATTTCGATGGTGCGGACCTGGTTTTCCATTTCTTTCATGTCTGTTTCATCATCCCATGGCTTGACGTCAAGGATTATAGATGATTTAGCAATGAGAGCAGGTTTCTTCGATTTTTTGTCAGCATAAGCTTTCAGACGTTCTTCACGAATTCTGGCGGCCTCTGCATCCTCTTCCTCGTCGCCAGATCCGAATAGGTCGACGTCATCATCGTCGTCATCCTTAGCAGCGGGTGCAGGTGCTGTGGTGGTCTTTCCGCCAGCGGTCAGAGGGCTAGCTCCAGCCGCCCAGGATTTGCGCTCCTCGGCATTGTAGGAGGCAATCTGGTTGTACCACCGGAGAGCGTGGGGCAAGCTTGCGCTCGGGGCCTTGCCGACTTGTTGGAACACTTGGAAATCGGCTTGAGTTGGTGTGTAACCTGACACATAACTCTTATCAGCTAAGAAATTGTTAAGCTCATTAAGGCCTTGTGCGGTTTTAACGTCTCCTACTGCCATTTTTATTTAAAAAAGTAGTTTTTTTCGAGTAATAGCCACACGCCCCGGCACAATGTGAGCAAGAAGGAATGAAAAAGAAATCTGACATTGACATTGCCATGAAATTGACTTTCAAAGAACGAATGAATTGAACTAATTTGAACGG

Я только получаюжелаемый вывод для сотого идентификатора из моего ids.txt.Может ли кто-нибудь помочь мне, где мой сценарий не так.Я хотел бы получить все 100 последовательностей при запуске сценария.Спасибо

Я добавил ссылки на диски Google в файлы, с которыми я работаю: ids.txt

Source.fasta

Answer 1

лучший способ сделать это - использовать Python или Perl.Мне удалось создать скрипт для извлечения идентификаторов с помощью python следующим образом.

#script to extract sequences from a source file based on ids in another file
#the source is a fasta file with a header and a sequence that follows in one line
#the ids file contains one id per line
#both the id and source file should contain the character '>' at the beginning that siginifies an id

def main():

    #asks the user for the ids file 
    file1 = raw_input('ids file: ');
    #opens the ids file into the memory
    ids_file = open(file1, 'r');
    #asks the user for the fasta file
    file2 = raw_input('fasta file: ');
    #opens the fasta file into memory; you need your memory to be larger than the filesize, or python will hard crash
    fasta_file = open(file2, 'r');

    #ask the user for the file name of output file
    file3 = raw_input('enter the output filename: ');
    #opens output file with append option; append is must as you dont want to override the existing data
    output_file = open(file3, 'w');

    #split the ids into an array
    ids_lines = ids_file.read().splitlines()
    #split the fasta file into an array, the first element will be the id followed by the sequence
    fasta_lines = fasta_file.read().splitlines()

    #initializing loop counters
    i = 0;
    j = 0;

    #while loop to iterate over the length of the ids file as this is the limiter for the program
    while j<len(fasta_lines) and i<len(ids_lines):
            #if statement to match ids from both files and bring matching sequences
            if ids_lines[i] == fasta_lines[j]:
                #output statements including newline characters
                output_file.write(fasta_lines[j])
                output_file.write('\n')
                output_file.write(fasta_lines[j+1])
                output_file.write('\n')
                #increment i so that we go for the next id
                i=i+1;
                #deprecate j so we start all over for the new id
                j=0;
            else:
                #when there is no match check the id, we are skipping the sequence in the middle which is j+1
                j=j+2;

    ids_file.close()
    fasta_file.close()
    output_file.close()

main()`

Код не идеален, но работает для любого количества идентификаторов.Я проверил свои образцы, которые содержали 5000 идентификаторов в одном из них, и программа работала нормально.Если есть улучшения в коде, пожалуйста, сделайте это, я относительно новичок в программировании, поэтому код немного грубоват.

Answer 2

Это может работать для вас (GNU sed):

sed 's#.*#/^&/{s/ .*//;N;p}#' idFile | sed -nf - fastafile

Преобразовать idFile в скрипт sed и запустить его для fastaFile.

Answer 3

$ awk 'NR==FNR{a[$1];next} $1 in a{c=2} c&&c--' ids.txt source.fasta
>TRINITY_DN80_c0_g1_i1 len=723 path=[700:0-350 1417:351-368 1045:369-722] [-1, 700, 1417, 1045, -2]
CGTGGATAACACATAAGTCACTGTAATTTAAAAACTGTAGGACTTAGATCTCCTTTCTATATTTTTCTGATAACATATGGAACCCTGCCGATCATCCGATTTGTAATATACTTAACTGCTGGATAACTAGCCAAAAGTCATCAGGTTATTATATTCAATAAAATGTAACTTGCCGTAAGTAACAGAGGTCATATGTTCCTGTTCGTCACTCTGTAGTTACAAATTATGACACGTGTGCGCTG

Выше было выполнено с использованием вашего опубликованного source.fasta и этого ids.txt:

$ cat ids.txt
>TRINITY_DN14840_c10_g1_i1
>TRINITY_DN80_c0_g1_i1

Answer 4

Первая группа всех идентификаторов, как одно выражение, разделенных |как это

cat ids.txt | tr '\n' '|' | awk "{print "\"" $0 "\""}'

Удалить последний |символ из выражения.

Теперь вы можете выполнить grep, используя выходные данные, полученные из предыдущей команды, например:

egrep -E ">TRINITY_DN14840_c10_g1_i1|>TRINITY_DN8506_c0_g1_i1|>TRINITY_DN12276_c0_g2_i1|>TRINITY_DN15434_c5_g3_i1|>TRINITY_DN9323_c8_g3_i5|>TRINITY_DN11957_c1_g7_i1|>TRINITY_DN15373_c1_g1_i1|>TRINITY_DN22913_c0_g1_i1|>TRINITY_DN13029_c4_g5_i1" source.fasta

Это будет печатать только совпадающие строки

Редактирование согласноtripleee comments

Использование следующего - правильная печать вывода. Предположим, что ID и последовательность находятся на разных линиях

tr '\n' '|' <ids.txt | sed 's/|$//' | grep -A 1 -E -f - source.fasta

Answer 5

Повторное зацикливание большого файла неэффективно;вы действительно хотите избежать запуска grep (или sed или awk) более одного раза, если вы можете избежать этого.Вообще говоря, sed и Awk часто легко позволяют вам указать действия для отдельных строк в файле, а затем запустить скрипт над файлом только один раз.

Для этой конкретной проблемы используется стандартная идиома Awk.с NR==FNR пригодится.Это механизм, который позволяет считывать несколько ключей в память (конкретно, когда NR==FNR это означает, что вы обрабатываете первый входной файл, потому что общий номер строки ввода равен номеру строки в этом файле), а затемпроверьте, присутствуют ли они в последующих входных файлах.

Напомним, что Awk читает по одной строке за раз и выполняет все действия, условия которых совпадают.Условия представляют собой простое логическое значение, а действия представляют собой набор команд Awk в паре фигурных скобок.

awk 'NR == FNR { s[$0]; next }
    # If we fall through to here, we have finished processing the first file.
    # If we see a wedge and p is 1, reset it -- this is a new sequence
    /^>/ && p { p = 0 }
    # If the prefix of this line is in s, we have found a sequence we want.
    ($1$2 in s) || ($1 in s) || ((substr($1, 1, 1) " " substr($1, 2)) in s) {
        if ($1 ~ /^>./) { print $1 } else { print $1 $2 }; p = 1; next }
    # If p is true, we want to print this line
    p' ids.txt source.fasta >out.fasta

Поэтому, когда мы читаем ids.txt, условие NR==FNR выполняется, и поэтомумы просто храним каждую строку в массиве s.next заставляет пропустить оставшуюся часть сценария Awk для этой строки.

При последующих чтениях, когда NR!=FNR, мы используем переменную p для управления тем, что печатать.Когда мы видим новую последовательность, мы устанавливаем p в 0 (если это было 1 из предыдущей итерации).Затем, когда мы видим новую последовательность, мы проверяем, находится ли она в s, и если да, мы устанавливаем p в единицу.Последняя строка просто печатает строку, если p не пусто или равно нулю.(Пустое действие является сокращением для действия { print }.)

Немного сложное условие, чтобы проверить, находится ли $1 в s, может быть слишком сложным - я ввел некоторые нормализации, чтобы убедиться,что пробел между > и идентификатором последовательности допускается, независимо от того, был ли он в ids.txt или нет.Вероятно, это можно упростить, если ваши файлы последовательно отформатированы.

Answer 6

Только с GNU grep и sed:

grep -A 1 -w -F -f ids.txt source.fasta | sed 's/ .*//'

---

См .: man grep