我想用linux shell统计人类参考基因组中每个染色体的长度和所有染色体的总长度,编写bash脚本如下:
#!/bin/bash
# 函数用于计算每个染色体的总长度并保存到指定的输出文件中
# 使用方法:calculate_chromosome_lengths input_reference_file output_file
function calculate_chromosome_lengths {
local reference_file="$1" # 输入的FASTA文件路径
local output_file="$2" # 输出文件路径
# 声明关联数组用于存储每个染色体的长度
declare -gA chromosome_lengths
echo "Chromosome Chromosome_Length" > a.txt
# 使用seqkit将FASTA文件转换为TAB格式,并提取染色体名称和序列信息
seqkit fx2tab "$reference_file" | awk '{print $1"\t"$NF}' | grep -v chrKI | grep -v chrGL |
while IFS=$'\t' read -r chromosome sequence; do
# 统计染色体的长度
length=${#sequence}
chromosome_lengths["$chromosome"]=$length
echo "Debug: "$chromosome": ${chromosome_lengths["$chromosome"]}"
echo ""$chromosome" ${chromosome_lengths["$chromosome"]}" >> a.txt
done
echo "Debug: Length of chr1: ${chromosome_lengths[chr1]}"
# 将结果输出到指定的输出文件中
echo "Chromosome Chromosome_Length" > "$output_file"
for chromosome in "${!chromosome_lengths[@]}"; do
echo "$chromosome ${chromosome_lengths[$chromosome]}" >> "$output_file"
done
# 统计所有染色体的总长度,并将总长度添加在输出文件的最后一行
total_length=0
for length in "${chromosome_lengths[@]}"; do
total_length=$((total_length + length))
done
echo "Total $total_length" >> "$output_file"
echo "染色体长度统计已完成,并保存到 $output_file 文件中。"
}
# 使用示例:将FASTA文件作为第一个参数,输出文件作为第二个参数
if [[ $# -lt 2 ]]; then
echo "请提供输入的FASTA文件和输出文件路径。"
echo "使用方法:$0 input_reference_file output_file"
else
calculate_chromosome_lengths "$1" "$2"
fi
我目前遇到的问题是echo "Debug: Length of chr1: ${chromosome_lengths[chr1]}" 这一步之前,数组chromosome_lengths输出正常,但是chromosome_lengths并没有传递到下一步,echo "Debug: Length of chr1: ${chromosome_lengths[chr1]}"的结果为: "Debug: Length of chr1:"。输出结果为空,请问为什么?