Вопросы Теги

Извлечение последовательностей из файла fasta

У меня есть файл fasta (не в правильном формате), который содержит сотни тысяч последовательностей ДНК разной длины, как это:

>NODE_213384_length_62_cov_8686_ID_2134025ATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAACGGAATCGAATGG>NODE_213385_length_62_cov_7933_ID_2134027ATCATCATCGAATGGAATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAAC>NODE_213386_length_62_cov_7184_ID_2134029AATGATTATTCCATTCGAGTCCATTCGATGATTCCATTCGATTCCATTCGATGATGATTGCA>NODE_213387_length_62_cov_8639_ID_2134031CAGAGCAGACTTGAAACACTCTTTTTGTGGAATTTGCAAGTGGAGATTTCAGCCGCTTTGAG>NODE_213388_length_62_cov_6833_ID_2134033AGACTTGAAACACTCTTTTTGTGGAATTTGCAAGTGGAGATTTCAGCCGCTTTGAGGTCAAT

Я хотел бы использовать простую команду Linux для извлечения последовательностей длиннее 1000bp и вывода в правильный формат fasta, как это:

>NODE_213384_length_62_cov_8686_ID_2134025
ATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAACGGAATCGAATGG
>NODE_213385_length_62_cov_7933_ID_2134027
ATCATCATCGAATGGAATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAAC
>NODE_213386_length_62_cov_7184_ID_2134029
AATGATTATTCCATTCGAGTCCATTCGATGATTCCATTCGATTCCATTCGATGATGATTGCA
>NODE_213387_length_62_cov_8639_ID_2134031
CAGAGCAGACTTGAAACACTCTTTTTGTGGAATTTGCAAGTGGAGATTTCAGCCGCTTTGAG

Буду благодарен любому, кто сможет помочь с этим.

-2
15.03.2016, 22:44
3 ответа

Вы можете сделать это в два этапа; сначала разделите на целевой формат, затем распечатайте пары строк, в которых последовательность ДНК достаточно длинная. Например, если предположить, что 1000 пар оснований относится к длине символа последовательности ДНК (а не к значению после слова «длина»), это может быть следующее 

cat inputfile | sed -e 's/>/\n>/g' -e 's/\(ID_[0-9]*\)/\n\1/g' |\
awk -F_ '$1=="NODE" { n=$0; next; } length($0)>1000 { print n; print ">" $0;}' \
> outputfile
0
28.01.2020, 05:18

Простая команда Linux" может выглядеть примерно так:

sed 's/\(>[^ACGT]*_[0-9]\+\)\([ACGT]\+\)/\1\n\2\n/g' yourdnafile |egrep -B1 '^[ACGT]{1000}'

Часть sed разбивает на 2 строки на набор, а grep показывает совпадение более 1000 и строки, предшествующей этой (-B1).

или это может быть еще проще:

sed 's/\(>[^>ACGT]*\)\([ACGT]\{1000\}[ACGT]*\)/\1\n\2\n/g;s/>[^>ACGT\n]*[ACGT]\+//g' yourdnafile
0
28.01.2020, 05:18

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

head -20 inputfile | 
tr '>' '\n'  > stage1
perl -ne 'print ">$1 $2\n" if /^(.*?)([ACGTU]+)$/ && length($2)>1000' < stage1 > output

Когда вы убедитесь, что все работает с первыми 20 строками, сделайте это по-настоящему в одном конвейере: 

tr '>' '\n' inputfile | 
perl -ne 'print ">$1 $2\n" if /^(.*?)([ACGTU]+)$/ && length($2)>1000' > output

Мой perl-скрипт, возможно, не так эффективен, как другие, но должен справиться с задачей. Я написал его для наглядности. Выведите метку, за которой следует пробел, затем соответствующие пары оснований и новую строку, если и только если есть такая строка, содержащая более 1000 пар оснований.

0
28.01.2020, 05:18

Теги

Похожие вопросы