如何写shell脚本统计fasta序列的长度和gc含量

可以使用以下shell脚本来统计fasta序列的长度和GC含量：

#!/bin/bash

# 读取fasta文件
fasta_file=$1
fasta=$(cat $fasta_file)

# 统计序列长度
len=$(echo "$fasta" | grep -v "^>" | tr -d "\n" | wc -m)
echo "Sequence length: $len"

# 统计GC含量
gc=$(echo "$fasta" | grep -v "^>" | tr -d "\n" | grep -o "[GCgc]" | wc -l)
gc_perc=$(echo "scale=2;$gc/$len*100" | bc)
echo "GC content: $gc_perc%"

使用该脚本时，需要将fasta文件作为第一个参数传递给脚本，例如：

bash fasta_stats.sh input.fasta

该脚本首先读取fasta文件，然后使用grep和tr命令来统计序列长度和GC含量。最后，使用echo命令打印结果。

相关问题

提取fasta序列中的指定位置的长度序列的脚本

以下是一个 Python 脚本，可用于从 FASTA 格式的文件中提取指定位置的长度序列：

def read_fasta_file(filename):
    """读取 FASTA 格式的文件并返回一个字典，其中键为序列 ID，值为序列字符串"""
    sequences = {}
    with open(filename, 'r') as f:
        current_id = None
        current_sequence = ''
        for line in f:
            line = line.strip()
            if line.startswith('>'):
                if current_id is not None:
                    sequences[current_id] = current_sequence
                current_id = line[1:]
                current_sequence = ''
            else:
                current_sequence += line
        if current_id is not None:
            sequences[current_id] = current_sequence
    return sequences

def extract_sequence(sequence, start, end):
    """从序列中提取指定位置的长度序列"""
    return sequence[start-1:end]

if __name__ == '__main__':
    # 读取 FASTA 文件
    sequences = read_fasta_file('sequences.fasta')

    # 提取指定位置的长度序列
    sequence_id = 'sequence1'
    start = 10
    end = 20
    extracted_sequence = extract_sequence(sequences[sequence_id], start, end)

    # 输出结果
    print(f'>{sequence_id}_{start}_{end}')
    print(extracted_sequence)

你只需要将上述脚本中的 `sequences.fasta` 文件名、`sequence1` 序列 ID、`start` 和 `end` 位置替换为你自己的需要即可。

shell编程读取一个fasta序列，输出小于氨基酸序列小于150的序列

你可以使用下面的shell脚本来实现：

#!/bin/bash

while read line
do
    if [[ ${line:0:1} == '>' ]]
    then
        if [[ ! -z $seq && $seqlen -le 150 ]]
        then
            echo -e "$header\n$seq"
        fi
        header=$line
        seq=''
        seqlen=0
    else
        seq=$seq$line
        seqlen=$((seqlen+${#line}))
    fi
done < your_sequence.fasta

if [[ $seqlen -le 150 ]]
then
    echo -e "$header\n$seq"
fi

这个脚本会逐行读取 `your_sequence.fasta` 中的内容。对于每一行，如果它以 `>` 开头，则表示前面的序列已经读完，需要对前一个序列进行处理。如果前一个序列的长度小于等于 150，则输出该序列。如果当前行不是以 `>` 开头，则说明当前行是序列的一部分，将其加入 `seq` 变量中，并更新序列的长度 `seqlen`。最后，当文件读取完毕时，如果最后一个序列的长度小于等于 150，则输出该序列。

需要注意的是，这个脚本假设每个序列都以一个以 `>` 开头的行开始，并且序列中没有空行。如果您的数据有所不同，您可能需要对脚本进行适当的修改。