다음과 같은 파일이 있습니다.
CDS join(36..56,37..67)
CDS 36..183
CDS 457..565
CDS join(505..519,521..596)
CDS join(577..591,725..770)
CDS join(516..591,725..899)
CDS 508..556
CDS 571..841
CDS complement(619..788)
CDS 843..863
파일에 있는 특정 수의 뉴클레오티드 범위를 인쇄하고 싶습니다(다른 파일 "sequence.fasta"에서 시퀀스를 읽습니다). 예를 들어, Sequence.fasta 파일은 다음과 같습니다.
>gi1234 HIVgenome|NC_909999.1
AACTGCGTGTGTGTCCACACAACACTGGGGGACACACAACAACAACACTGGGGGACACACTGGGACAACACTGGGGGACAGGACACTGTACAACACTGGGTGTGTCGGGACAGTACACATGTTGGGGGGGTGTGTCGGACAACACTGGGGGACATGTGTGTACAACACTGGGGGACAGTGACGACGACAACACTGGGGGACACGAGCGTTGTGAGCAGGTGACAACACTGGGGGACAGTGTTTTTACAACACTGGGGGACATTTTTGAGCAGCGACGCAGCGTTGTGGGGTGTGTCGGAAGGTGTGTCGTGTGTCGTGTGTC
출력 p는 다음과 같아야 합니다.
36 - 56 ACAACAACAACACTGGGGGAC
37 - 67 CAACAACAACACTGGGGGACAACACTGGGAC
& 곧...
~까지
843 - 863 GTGT....
쉘 스크립트를 통해 이를 수행하는 가장 쉬운 방법은 무엇입니까?
답변1
이 문제를 해결하려면 이 포럼에서 제공하는 것보다 더 많은 프로그래밍 노력이 필요합니다(저는 이런 종류의 프로그래밍을 직업으로 하고 있습니다).
이것DDBJ/ENA/GenBank 파일 형식(질문의 첫 번째 파일)은 복잡하며 CDS(게놈 서열의 코딩 부분)가 단순하거나 연결될 뿐만 아니라 상호보완적이며 이들의 조합이 될 수 있습니다. 또한 위치 좌표는 다음을 가질 수 있습니다.수정자일반적인 해결책을 위해서는 이 문제를 처리해야 합니다.
지역 생물정보학자(또는 프로그래머)에게 문의하거나 StackExchange와 같은 생물정보학 포럼에 문의하는 것이 더 나을 것입니다.생물정보학 웹사이트. 이런 종류의 작업을 수행하기 위한 기존 도구를 알려주거나, 생물정보학자를 알고 있으면 더 효율적일 수 있는 이상한 BioPerl/BioPython 스크립트를 제공할 것입니다 ;-)
하나가능한경로를 사용 중입니다.GenBank 특징 추출기, 그러나 온라인에서 사용하는 것은 소규모 데이터 세트 이외의 경우에는 최선의 선택이 아닐 수 있습니다.
답변2
이전 답변에도 불구하고 fasta 파일에서 특정 하위 시퀀스를 추출하는 하위 문제를 조사했습니다. 솔루션은 두 부분으로 나뉩니다.
sh
명령줄 구문 분석 및 호출을 수행하는 쉘 스크립트 ...awk
fasta 파일을 구문 분석하는 스크립트입니다.
보여주기 때문에 여기에 게시하기로 결정했습니다.
- 쉘 스크립트에서 옵션의 명령줄 구문 분석을 수행하는 방법.
- 단순한 “한 줄”
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;
}
답변3
비슷한 문제를 해결하고 있지만 과거에 답변(WSL, Mac) 중 하나에 제공된 스크립트를 사용하려고 할 때 약간의 문제가 있었습니다...
그래서 저는 계속 찾아보았고 매우 간단한 해결책을 찾았습니다. 그것이 바로 키트 getfasta
에 포함되어 있는 것이었습니다 bedtools
!
https://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html