如何在R语言中基于读取数剔除低读取数Fastq文件？

方法1：用系统自带工具实现（无需额外安装） #

Fastq格式的每条读取固定占4行，因此读取数等于文件总行数除以4。可以通过命令行循环遍历文件，计算行数后判断是否达标：

# 遍历当前目录下所有fastq/fastq.gz文件（可根据实际后缀调整）
for file in *.fastq *.fastq.gz; do
    # 区分压缩与非压缩文件，计算总行数
    if [[ $file == *.gz ]]; then
        line_count=$(zcat "$file" | wc -l)
    else
        line_count=$(cat "$file" | wc -l)
    fi
    # 计算读取数
    read_count=$((line_count / 4))
    # 筛选读取数≥10的文件，可选择写入列表或直接移动到目标目录
    if [[ $read_count -ge 10 ]]; then
        echo "$file" >> valid_samples.txt
        # 可选操作：将符合条件的文件移动到指定文件夹
        # mv "$file" ./valid_samples/
    fi
done

这个方案依赖bash、wc等系统默认工具，适合环境受限的场景。

方法2：用seqkit工具（高效直观） #

如果允许安装第三方工具，seqkit是专门处理序列文件的轻量工具，能直接输出Fastq的读取数，操作更简洁：

# 生成所有文件的统计信息，筛选读取数≥10的文件
seqkit stats *.fastq *.fastq.gz | awk 'NR>1 && $4>=10 {print $1}' >> valid_samples.txt

解释：seqkit stats会输出文件路径、总序列数（即读取数）等信息，awk跳过表头（NR>1），判断第4列（序列数）≥10的行并提取文件名。

补充提示 #

• 若处理双端测序的配对文件，需注意同步保留或剔除配对的两个文件：可先筛选单端有效文件，再根据文件名规则（如_R1/_R2）匹配对应的配对文件。
• 对于格式不规范的Fastq文件（含空行等），可先用seqkit seq -c "$file"验证行数是否为4的倍数，避免误判读取数。

内容的提问来源于stack exchange，提问作者Sara Nicholson