☰

목록에서 fasta 항목을 추출하려면 읽는 동안 사용하세요.

tag-icon • tag-icon text-processing sed bioinformatics
목록에서 fasta 항목을 추출하려면 읽는 동안 사용하세요.

나는 각각 약 14,000개의 "항목"을 포함하는 28개의 파일을 가지고 있습니다. 단일 항목은 헤더(>string으로 표시됨), 개행 문자 및 일련의 문자열로 구성됩니다. 각 항목에는 가변 길이 시퀀스/문자열이 있습니다. 28개 파일은 모두 동일한 항목 헤더를 갖고 있지만 각 항목의 순서는 다양합니다.

예를 들어 CR1_ref.fasta 파일은 다음과 같습니다.

>FBgn0080937
ATGGATAAAAGGCTCAGCGATAGTCCCGGAGATTGTCGCGTAACCAGATCCAGCATGACGCCCACCCTCCGCTTGGAGCACAGTCCCCGGCGGCAACAACAGCAACAACA
>FBgn0076379
ATGCTGCGCACCCTTTTCGCCGTGCGTGGTCAGTGCCAGCAGCTGCTGAGGAGAACATTCACCCCCCATTGCAGTGGCCAACGA
>FBgn0070974
ATGCAGACGCGTCCGAGCAGTGAACCGCAGCGCGCCAAGGAGCAACTCCTGCGGGAGCTGCCGCCGCAGAAATGCTCCAGCGCCACGCTGGCCAAGAAGGTGCTGTCGCAGAGCCCGCCGGCAGCCCCGCCGCCCACACCGGCCACAATTGTGCCGCTCACTGCGGTGCCCGTCATCCAGCTGACGCCTCCGTCGCACTCCGGCGACACGCCGCAAAAGCCAGCACCTCCGGCGCCGCCGCCGCC

전반적인 목표는 약 14,000개의 새 파일을 생성하는 것입니다. 이러한 각 파일은 28개 파일 전체의 특정 ID/헤더와 연결된 항목입니다.

단일 파일에서 단일 항목을 추출하려면 다음 명령을 사용할 수 있습니다.

sed -n '/^>FBgn0080937$/{p;n;p;}' CR1_ref.fasta

28개 파일 모두(각 파일은 ref.fasta로 끝남)에서 이 항목을 추출하려면 다음을 수행할 수 있습니다.

for i in *ref.fasta; do sed -n '/^>FBgn0080937$/{p;n;p}' $i; done > FBgn0080937.fasta

나는 gene.txt라는 항목의 헤더에 각각 해당하는 14,000줄로 구성된 별도의 텍스트 파일을 가지고 있습니다. 파일의 처음 몇 줄은 다음과 같습니다.

FBgn0080937
FBgn0076379
FBgn0070974
FBgn0081668
FBgn0076576
FBgn0076572
FBgn0079684
FBgn0070907
FBgn0080226
FBgn0072746

파일을 읽고 각 헤더 ID에 대해 새 텍스트 파일을 만들고 싶습니다. 아래 $F는 특정 헤더(FBgn*)에 대한 항목을 추출하여 새 파일에 저장합니다. 나는 대체 명령을 사용하여 해당 시퀀스의 원본인 while ref.fasta 파일을 기반으로 시퀀스의 이름을 바꾸고 있습니다.

while read -r line;
do F=$line
for i in *ref.fasta
do sed -n "/^>$F$/{s/FB.*/$i/;p;n;p;}" $i > $line.fasta
done
done < "gene.txt"

현재 스크립트는 14,000개의 파일을 생성하지만 각 파일에는 하나의 시퀀스만 있습니다.

>Z9_ref.fasta
ATGCAGACGCGTCCGAGCAGTGAACCGCAGCGCGCCAAGGAGCAAC

각 *ref.fasta 파일에는 시퀀스당 28개의 시퀀스가 ​​있을 것으로 예상됩니다. sed 명령이 마지막 항목을 출력하고 있습니다. 예상 출력은 다음과 같습니다.

    >CR1_ref.fasta
    ATGCAGACGCGTCCGAGCAGTGAACC
    >FH2_ref.fasta
    AGCAGTGAACCGCAGCGCGCCAAGGAGCAAC
    >MSH10_ref.fasta
    CGCGTCCGAGCAGTGAACCGCAGCGCGCCAAGGAGCAAC
    >Z9_ref.fasta
    ATGCAGACGCGTCCGAGCAGTGAACCGCAGCGCGCCAAGGAGCAAC

답변1

쉘은 이러한 유형의 구문 분석에 실제로 적합하지 않습니다. 코드에서 전체 파일을 한 번 읽는 것을 볼 수 있습니다.각파일에서 읽은 유전자 이름입니다 gene.txt.

아래의 단일 명령은 awk동일한 작업을 더 빠르게 수행합니다.

awk -F '>' '
    FNR == NR           { genes[$1]; next }
    /^>/ && $2 in genes { if (out != "") close(out);
                          out = $2 ".fa"
                          split(FILENAME, a, "_")
                          $0 = ">" a[1] "_" $2 }
    out != ""           { print >>out }' genes.txt *_ref.fasta

먼저 genes.txt파일을 읽고 genes유전자 이름을 키로 사용하여 해당 파일에서 호출되는 연관 배열을 만듭니다.

Fasta 파일에 도달하면(코드에서는 이러한 파일이 모두 다음과 같이 호출된다고 가정합니다 XXX_ref.fasta), Fasta 헤더를 읽고 헤더의 유전자가 목록의 키일 genes때 유전자 이름에서 출력을 생성합니다. filename을 genename.fa밑줄 앞에 현재 파일 이름 부분을 포함하도록 헤더를 다시 작성합니다.

원래 헤더가 XXX_ref.fasta다음과 같은 경우

>genename

그러면 이것은 다음으로 변환됩니다.

>XXX_genename

스크립트의 마지막 부분은 awk모든 행을 적절한 출력 파일로 보냅니다.

제공한 데이터로 테스트하면 다음 세 가지 파일이 생성됩니다.

$ ls *.fa
FBgn0070974.fa FBgn0076379.fa FBgn0080937.fa

$ cat FBgn0076379.fa
>CR1_FBgn0076379
ATGCTGCGCACCCTTTTCGCCGTGCGTGGTCAGTGCCAGCAGCTGCTGAGGAGAACATTCACCCCCCATTGCAGTGGCCAACGA

관련 정보