☰

fasta 파일에서 시퀀스 추출

tag-icon • tag-icon text-processing
fasta 파일에서 시퀀스 추출

다음과 같이 다양한 길이의 수십만 개의 DNA 서열을 포함하는 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

도움을 주실 수 있는 분에게 감사드립니다.

답변1

줄이 긴 것 같네요. sed 및 awk가 이 작업을 처리하기 위해 사용 가능한 메모리를 사용하므로 문제가 발생할 수 있습니다(물론 파일 크기/줄 길이에 따라 다름). 따라서 2단계 접근 방식이 채택되지만 메모리 제한으로 인해 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개 이상의 염기쌍을 포함하는 줄이 있는 경우에만 라벨을 인쇄하고 그 뒤에 공백, 관련 염기쌍 및 개행 문자를 인쇄합니다.

답변2

"간단한 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

답변3

먼저 대상 형식으로 분할한 다음 DNA 서열의 충분히 긴 라인 쌍을 인쇄하는 두 단계로 이 작업을 수행할 수 있습니다. 예를 들어, 1000bpDNA 서열 문자 길이("length"라는 단어 뒤의 값이 아닌)를 참조한다고 가정하면 다음과 같을 수 있습니다.

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

관련 정보