fasta 파일에서 하위 집합 추출

Question 1

모든 시퀀스가 한 줄에 맞는 간단한 예의 경우 다음을 사용할 수 있습니다 grep( grep해당 옵션을 지원하지 않는 경우 --no-group-separator출력을 전달 grep -v -- '--').

$ grep -wEA1 --no-group-separator 'chr1|chr2|chr21|chrX' file.fa 
>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

여러 행의 시퀀스가 있다고 가정하면(염색체를 다루는 경우 그럴 가능성이 높음) 먼저 이를 행으로 연결해야 합니다. 이로 인해 상황이 상당히 복잡해졌습니다. awk한 줄을 사용할 수 있습니다 .

$ awk -vRS=">" 'BEGIN{t["chr1"]=1;t["chr2"]=1;t["chr21"]=1;t["chrX"]=1}
                {if($1 in t){printf ">%s",$0}}' file.fa
>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

위 스크립트는 레코드 구분 기호를 >( vRS=">")로 설정합니다. 이는 "선"이 >~및 에 의해 정의 됨을 의미합니다 \n. 그런 다음 블록은 BEGIN각 대상 시퀀스 ID가 키인 배열을 설정합니다. 나머지는 각 "라인"(시퀀스)을 확인하고 첫 번째 필드가 배열( )에 있으면 앞에 $i in t현재 "라인"( )을 인쇄합니다 .$0>

이런 종류의 작업을 자주 수행하면 배열 작성이 금방 짜증날 수 있습니다. 개인적으로 저는 FASTA를 tbl 형식( sequence_name<TAB>sequence한 줄에 하나의 시퀀스) 으로 변환하고 그 반대로 변환하기 위해 다음 두 스크립트(이전 연구실 파트너로부터 상속받았음)를 사용합니다 .

FastaToTbl

#!/usr/bin/awk -f
{
        if (substr($1,1,1)==">")
        if (NR>1)
                    printf "\n%s\t", substr($0,2,length($0)-1)
        else 
            printf "%s\t", substr($0,2,length($0)-1)
        else 
                printf "%s", $0
}END{printf "\n"}

TblToFasta

#! /usr/bin/awk -f
{
  sequence=$NF

  ls = length(sequence)
  is = 1
  fld  = 1
  while (fld < NF)
  {
     if (fld == 1){printf ">"}
     printf "%s " , $fld
     if (fld == NF-1){
        printf "\n"
      }
      fld = fld+1
  }
  while (is <= ls){
    printf "%s\n", substr(sequence,is,60)
    is=is+60
  }
}

파일에 저장 $PATH하고 실행 가능하게 만들면 grep대상 시퀀스에 간단히 사용할 수 있습니다(위와 달리 여러 줄 시퀀스에 대해 작동함).

$ FastaToTbl file.fa | grep -wE 'chr1|chr2|chr21|chrX' | TblToFasta
>chr1 
ACGGTGTAGTCG
>chr2 
ACGTGTATAGCT
>chr21 
ACGTTGATGAAA
>chrX 
GTACGGGGGTGG

grep검색 대상 파일을 전달할 수 있으므로 확장하기가 더 쉽습니다.

$ cat ids.txt 
chr1
chr2
chr21
chrX

$ FastaToTbl file.fa | grep -wFf ids.txt | TblToFasta
>chr1 
ACGGTGTAGTCG
>chr2 
ACGTGTATAGCT
>chr21 
ACGTTGATGAAA
>chrX 
GTACGGGGGTGG

마지막으로 대규모 시퀀스를 처리하는 경우 이를 수행할 수 있는 다양한 도구 중 하나를 사용하는 것을 고려할 수 있습니다. 장기적으로 보면 더 빠르고 효율적일 것입니다. 예를 fastafetch들어 ,exonerate모음곡. 대부분의 Linux 배포판의 리포지토리에서 사용할 수 있습니다. sudo apt-get install exonerate예를 들어 Debian 기반 시스템에서는 . 설치 후 다음을 수행할 수 있습니다.

## Create the index    
fastaindex -f file.fa -i file.in
## Loop and retrieve your sequences
for seq in chr1 chr2 chr21 chrX; do
     fastafetch -f file.fa -i file.in -q "$seq" 
done 
>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

또는 내 것을 사용할 수도 있습니다.retrieveseqs.pl, 여기에는 다른 멋진 기능도 있습니다.

$ retrieveseqs.pl -h

    retrieveseqs.pl will take one or more lists of ids and extract their sequences from 
    multi FASTA file

    USAGE : retrieveseqs.pl [-viofsn] <FASTA sequence file> <desired IDs, one per line>

    -v : verbose output, print a progress indicator (a "." for every 1000 sequences processed)
    -V : as above but a "!" for every desired sequence found.
    -f : fast, takes first characters of name "(/^([^\s]*)/)" given until the first space as the search string
         make SURE that those chars are UNIQUE.
    -i : use when the ids in the id file are EXACTLY identical
         to those in the FASTA file
    -h : Show this help and exit.
    -o : will create one fasta file for each of the id files
    -s : will create one fasta file per id
    -n : means that the last arguments (after the sequence file)
         passed are a QUOTED list of the names desired.
    -u : assumes uniprot format ids (separated by |)

귀하의 경우에는 다음을 수행합니다.

$ retrieveseqs.pl -fn file.fa "chr1 chr2 chr21 chrX"
[7 (4/4 found]
>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

이 내용은 제가 작업을 위해 작성한 것이며 제대로 문서화되거나 지원되지 않습니다. 그럼에도 불구하고 나는 몇 년 동안 그것을 행복하게 사용하고 있습니다.

Answer