递归计算多个 fastq 文件中的序列数

递归计算多个 fastq 文件中的序列数

我有许多以 fastq.gz 结尾的 fastq 文件。

rep1.fastq.gz
rep2.fastq.gz
rep3.fastq.gz
rep4.fastq.gz
.....

我期望我的输出为

rep1.fastq.gz 23516782
rep2.fastq.gz 45126780
rep3.fastq.gz 67543908
rep4.fastq.gz 76425368

其中第 1 行显示我的每个输入文件,第 2 行显示每个文件中的序列数。

为了实现这一点,我编写了一个小的 bash 脚本来计算每个文件中的序列数,并将每个文件后面写入的数字作为输出

for sample in *.fastq.gz;do echo -en $sample "\t";(zcat $sample|wc -l)/4|bc ;done

我收到错误:-bash:意外标记“/4”附近出现语法错误

答案1

忽略每个序列 4 行的假设中可能出现的所有错误...上面显示的命令应使用以下格式

for file in *.fastq.gz; do echo -en $file "\t";echo "$(zcat $file| wc -l)"/4 |bc;done

答案2

你正在运行这个:

(zcat $sample|wc -l)/4|bc

/4shell 无法理解其中唯一的内容,这就是它失败的原因。我怀疑你想要的是得到结果命令的值zcat $sample|wc -l,然后打印该值并将/4其传递给bc.如果是这样,您$()不仅需要而且()还需要引用它:

echo "$(zcat $sample|wc -l)/4" | bc

所以这意味着:

for sample in *.fastq.gz; do 
    echo -en $sample "\t"; echo "$(zcat $sample|wc -l)/4" | bc 
done

或者,更便携一点:

for sample in *.fastq.gz; do 
    printf '%s\t%s\n' "$sample" "$(echo "$(zcat "$sample" | wc -l)/4" | bc )"
done

或者,您可以在以下位置完成整个操作awk

for sample in *.fastq.gz; do 
    printf '%s\t' "$sample"
    zcat "$sample" | awk '!(NR % 4){k++}END{print k}'
done

但是,请注意,fastq 格式的定义中没有任何内容表明文件每个序列只有 4 行。如果你很了解你的数据,你可以使用这种方法,但如果你需要处理任意 fastq 文件,你不能假设只有 4 行,最好使用专用工具。

您可能会发现这个问答很有趣:快速计算 fastq 文件中读取次数和碱基数量的方法?

还有FASTQ文件格式规范澄清了您不能假设每个条目只有 4 行。也就是说,根据我过去 7 年在临床环境中处理人类 NGS 数据的丰富经验,我见过的每个文件每个样本只有 4 行。但我不处理长时间读取的数据,而且格式本身允许读取更多数据,因此这是需要考虑的事情。

相关内容