특정 숫자의 fasta 서열을 절단하고 ORF를 생성하는 방법

Question 1

이 질문에는 이 포럼에서 제공할 수 있는 것보다 더 큰 프로그래밍 노력이 필요합니다(저는 이런 종류의 프로그래밍을 생계로 하고 있습니다).

그만큼DDBJ/ENA/GenBank 파일 형식(문제의 첫 번째 파일)은 복잡하며 CDS(게놈 서열의 코딩 부분)가 단순하거나 결합된 것이 아니라 보완되고 조합될 수 있습니다. 또한 위치 좌표는 다음을 가질 수 있습니다.수정자일반적인 솔루션의 경우 이를 처리해야 합니다.

지역 생물정보학자(또는 프로그래머)에게 문의하거나 StackExchange와 같은 생물정보학 포럼에 문의하는 것이 더 나을 것입니다.생물정보학 사이트. 이런 종류의 작업을 수행하기 위한 기존 도구를 알려주거나, 생물정보학자를 알고 있으면 아마도 더 자주 작동할 기발한 BioPerl/BioPython 스크립트를 제공할 것입니다 ;-)

하나가능한경로는GenBank 특징 추출기, 그러나 온라인으로 사용하는 것은 소규모 데이터 세트 외에는 최선의 선택이 아닐 가능성이 높습니다.

Answer

이 질문에는 이 포럼에서 제공할 수 있는 것보다 더 큰 프로그래밍 노력이 필요합니다(저는 이런 종류의 프로그래밍을 생계로 하고 있습니다).

그만큼DDBJ/ENA/GenBank 파일 형식(문제의 첫 번째 파일)은 복잡하며 CDS(게놈 서열의 코딩 부분)가 단순하거나 결합된 것이 아니라 보완되고 조합될 수 있습니다. 또한 위치 좌표는 다음을 가질 수 있습니다.수정자일반적인 솔루션의 경우 이를 처리해야 합니다.

지역 생물정보학자(또는 프로그래머)에게 문의하거나 StackExchange와 같은 생물정보학 포럼에 문의하는 것이 더 나을 것입니다.생물정보학 사이트. 이런 종류의 작업을 수행하기 위한 기존 도구를 알려주거나, 생물정보학자를 알고 있으면 아마도 더 자주 작동할 기발한 BioPerl/BioPython 스크립트를 제공할 것입니다 ;-)

하나가능한경로는GenBank 특징 추출기, 그러나 온라인으로 사용하는 것은 소규모 데이터 세트 외에는 최선의 선택이 아닐 가능성이 높습니다.

Question 2

이전 답변을 제공했지만 fasta 파일에서 특정 하위 시퀀스를 추출하는 하위 문제를 살펴보았습니다. 해결책은 두 부분으로 구성됩니다.

sh명령줄 구문 분석을 수행하고 다음을 호출하는 쉘 스크립트입니다 .
awkfasta 파일을 구문 분석하는 스크립트입니다 .

보여주기 때문에 여기에 게시하기로 결정했습니다.

쉘 스크립트에서 옵션의 명령줄 구문 분석을 수행하는 방법.
단지 "한 줄짜리" awk가 아닌 스크립트를 작성하는 것이 가능하다는 것입니다 .awk

가정:

시퀀스 ID는 >헤더 줄 바로 뒤에 표시되며 그 뒤에 공백 문자가 옵니다.
시퀀스 데이터에는 공백이 없습니다.

스크립트는 이라고 합니다 extract.sh.

gi1234실행하려면 위치 36과 183(포함) 사이의 시퀀스 ID에 대한 시퀀스를 가져옵니다 .

$ sh extract.sh -i gi1234 -a 36 -b 183 <sequence.fa
ACAACAACAACACTGGGGGACACACTGGGACAACACTGGGGGACA
GGACACTGTACAACACTGGGTGTGTCGGGACAGTACACATGTTGGGGGGGTGTGTCGGACAACACTGGGGGACATGTGTG
TACAACACTGGGGGACAGTGACG

출력이 포맷되지 않았습니다. 이 경우에는 질문의 데이터를 가져와서 스크립트를 실행하기 전에 80번째 문자마다 줄 바꿈을 삽입했습니다.

쉘 스크립트( extract.sh):

#!/bin/sh

usage() {
    cat <<USAGE_END
Usage:
    sh $0 -i seq_id -a start -b end <sequence.fa
USAGE_END
}

while getopts "a:b:i:" opt; do
    case $opt in
        a) start_pos="$OPTARG"  ;;
        b) end_pos="$OPTARG"    ;;
        i) seq_id="$OPTARG"     ;;
        *) usage; exit 1        ;;
    esac
done

if [ "x$seq_id" = "x" ]; then
    echo "Missing sequence ID (-i seq_id)"
    usage
    exit 1
fi

if [ "x$start_pos" = "x" -o "x$end_pos" = "x" ]; then
    echo "Missing either start or end coordinates (-a, -b)"
    usage
    exit 1
fi

if [ ! -f "extract.awk" ]; then
    echo "Can not find the extract.awk AWK script!"
    exit 1
fi

awk -v id="$seq_id" -v a="$start_pos" -v b="$end_pos" -f extract.awk

스크립트 awk( extract.awk):

# state == 0:   not yet at wanted sequence entry in file
# state == 1:   found sequence entry, not yet at start position
# state == 2:   found start position, not yet at end position

# off:  offset in sequence read so far

BEGIN   { state = 0 }

state == 0 && $0 ~ "^>" id " " {
    # We found the sequence entry we're looking for.
    state = 1;
    off = 0;
    next;
}

state > 0 && /^>/  {
    # New sequence header found while parsing sequence, terminate.
    exit;
}

state == 1 {
    len = length($0);

    if (len + off < a) {
        # Start position is not on this line.
        off += len;
        next;
    }

    state = 2;

    if (len + off >= b) {
        # Whole sequence is found on this line.
        print(substr($0, a - off, b - a + 1))
        exit;
    }

    # Output partial start of sequence.
    print(substr($0, a - off , len - (a - off) + 1))

    off += len;

    next;
}

state == 2 {
    len = length($0);

    if (off > b) {
        # (we should not arrive here)
        exit;
    }

    if (len + off < b) {
        # End position is not on this line.
        print($0);
        off += len;
        next;
    }

    # We found the end position, output line until end position
    # and terminate.
    print(substr($0, 1, b - off));
    exit;
}

Answer