我有一个文件,如下所示:
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....
通过 shell 脚本实现这一点最简单的方法是什么?
答案1
这个问题需要比这个论坛提供的更大的编程工作量(我以这种编程为生)。
这DDBJ/ENA/GenBank 文件格式(问题中的第一个文件)很复杂,并且允许 CDS(基因组序列的编码部分)不仅是简单的或连接的,而且是互补的及其组合。此外,位置坐标可以有修饰语对于通用解决方案来说,需要处理这个问题。
您最好询问当地的生物信息学家(或程序员)或在生物信息学论坛,例如 StackExchange生物信息学网站。他们会向您指出用于执行此类操作的现有工具,或者,了解生物信息学家,为您提供一些古怪的 BioPerl/BioPython 脚本,这些脚本可能会更有效;-)
一可能的路线是使用GenBank 特征提取器,但是在线使用它很可能不是除了小数据集之外的任何东西的最佳选择。
答案2
尽管我给出了之前的答案,但我已经研究了从 fasta 文件中提取特定子序列的子问题。解决方案分为两部分:
- 一个
sh
执行命令行解析并调用...的 shell 脚本 awk
解析 fasta 文件的脚本。
我决定将其发布在这里,因为它表明
- 如何在 shell 脚本中对选项进行命令行解析。
- 可以编写一个
awk
脚本,而不仅仅是awk
“一句台词”。
假设:
- 序列 ID 将直接出现
>
在标题行的 后面,后跟一个空格字符。 - 序列数据中没有任何空格。
该脚本称为extract.sh
.
gi1234
要运行,请获取位置 36 和 183(含)之间的序列 ID :
$ sh extract.sh -i gi1234 -a 36 -b 183 <sequence.fa
ACAACAACAACACTGGGGGACACACTGGGACAACACTGGGGGACA
GGACACTGTACAACACTGGGTGTGTCGGGACAGTACACATGTTGGGGGGGTGTGTCGGACAACACTGGGGGACATGTGTG
TACAACACTGGGGGACAGTGACG
输出未格式化。在本例中,我获取了问题中的数据,并在运行脚本之前在每 80 个字符后插入换行符。
shell 脚本 ( 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