如何连接在不同车道上生成的RNA-seq文件

Question

我有很大的RNA-seq文件产生在不同的车道。我提取了几个文件名，如下所示。

MC9_FNEN_638A_S19_L008_R1_001.fastq.gz
MC9_FNEN_638A_S19_L008_R2_001.fastq.gz
MC9_FNEN_638A_S9_L001_R1_001.fastq.gz
MC9_FNEN_638A_S9_L001_R2_001.fastq.gz
MC9_FNEN_638A_S9_L002_R1_001.fastq.gz
MC9_FREN_638A_S9_L002_R2_001.fastq.gz
MC9_FREN_638A_S9_L006_R1_001.fastq.gz
MC9_FREN_638A_S9_L006_R2_001.fastq.gz
MC9_FREN_638A_S9_L008_R1_001.fastq.gz
MC9_FREN_638A_S9_L008_R2_001.fastq.gz
MC9_ZH_637A_S74_L001_R1_001.fastq.gz
MC9_ZH_637A_S74_L001_R2_001.fastq.gz
MC9_ZH_637A_S74_L003_R1_001.fastq.gz
MC9_ZH_637A_S74_L003_R2_001.fastq.gz
MC9_ZH_637A_S74_L007_R1_001.fastq.gz
MC9_ZH_637A_S74_L007_R2_001.fastq.gz
MC9_ZH_637A_S74_L008_R1_001.fastq.gz
MC9_ZH_637A_S74_L008_R2_001.fastq.gz
MC9_ZH_637A_S84_L008_R1_001.fastq.gz
MC9_ZH_637A_S84_L008_R2_001.fastq.gz
DR14_DCRP_479C_S50_L001_R1_001.fastq.gz
DR14_DCRP_479C_S50_L001_R2_001.fastq.gz
DR14_DCRP_479C_S50_L002_R1_001.fastq.gz
DR14_DCRP_479C_S50_L002_R2_001.fastq.gz
DR14_DCRP_479C_S50_L006_R1_001.fastq.gz
DR14_DCRP_479C_S50_L006_R2_001.fastq.gz
DR14_DCRP_479C_S50_L007_R1_001.fastq.gz
DR14_DCRP_479C_S50_L007_R2_001.fastq.gz
DR14_DCRP_479C_S50_L008_R1_001.fastq.gz
DR14_DCRP_479C_S50_L008_R2_001.fastq.gz

我想连接所有的序列在不同的车道产生的正读和反向读取。例如，前10行是来自同一动物和特定组织(MC9_FREN)的序列文件。我想连接在不同车道中生成的所有前读XXXXX_R1_001.fastq.gz，并将文件名MC9_FREN_R1.fastq.gz和所有反向读取XXXX_R2_001.fastq.gz到MC9_FREN_R2.fastq.gz。

cat MC9_FREN_638A_S19_L008_R1_001.fastq.gz MC9_FREN_638A_S9_L001_R1_001.fastq.gz  MC9_FREN_638A_S9_L002_R1_001.fastq.gz  MC9_FREN_638A_S9_L007_R1_001.fastq.gz MC9_FREN_638A_S9_L008_R1_001.fastq.gz > MC9_FREN_R1.fastq.gz
cat MC9_FREN_638A_S19_L008_R2_001.fastq.gz MC9_FREN_638A_S9_L001_R2_001.fastq.gz MC9_FREN_638A_S9_L002_R2_001.fastq.gz  MC9_FREN_638A_S9_L007_R2_001.fastq.gz MC9_FREN_638A_S9_L008_R2_001.fastq.gz  > MC9_FREN_R2.fastq.gz
cat MC9_ZH_637A_S74_L001_R1_001.fastq.gz MC9_ZH_637A_S74_L003_R1_001.fastq.gz  MC9_ZH_637A_S74_L007_R1_001.fastq.gz MC9_ZH_637A_S74_L008_R1_001.fastq.gz  MC9_ZH_637A_S84_L008_R1_001.fastq.gz > MC9_ZH_R1.gz
cat MC9_ZH_637A_S74_L001_R2_001.fastq.gz  MC9_ZH_637A_S74_L003_R2_001.fastq.gz MC9_ZH_637A_S74_L007_R2_001.fastq.gz MC9_ZH_637A_S74_L008_R2_001.fastq.gz MC9_ZH_637A_S84_L008_R2_001.fastq.gz > MC9_ZH_R2.gz
cat DR14_DCRP_479C_S50_L001_R1_001.fastq.gz DR14_DCRP_479C_S50_L002_R1_001.fastq.gz DR14_DCRP_479C_S50_L006_R1_001.fastq.gz DR14_DCRP_479C_S50_L007_R1_001.fastq.gz DR14_DCRP_479C_S50_L008_R1_001.fastq.gz  > DR14_DCRP_R1.gz   
cat DR14_DCRP_479C_S50_L001_R2_001.fastq.gz DR14_DCRP_479C_S50_L002_R2_001.fastq.gz  DR14_DCRP_479C_S50_L006_R2_001.fastq.gz DR14_DCRP_479C_S50_L007_R2_001.fastq.gz DR14_DCRP_479C_S50_L008_R2_001.fastq.gz  > DR14_DCRP_R1.gz

Answer 1

下面的循环为当前目录中的FastQ文件提供了唯一的文件名前缀。它依赖于这样一个事实:在我们想要的文件名前缀和文件名后面的R1或R2之间总是有四个下划线(D2)。

for name in *.fastq.gz; do
    printf '%s\n' "${name%_*_*_*_R[12]*}"
done | uniq

以下内容是等价的，但不使用循环(而不是删除文件名的最后一部分，而是保留文件名的第一部分)：

printf '%s\n' *.fastq.gz | sed 's/^\([^_]*_[^_]*\).*/\1/' | uniq

在给定的文件列表中，上述任何一个返回

DR14_DCRP
MC9_FNEN
MC9_FREN
MC9_ZH

然后，我们读取这些前缀并创建连接文件：

for name in *.fastq.gz; do
    printf '%s\n' "${name%_*_*_*_R[12]*}"
done | uniq |
while read prefix; do
    cat "$prefix"*R1*.fastq.gz >"${prefix}_R1.fastq.gz"
    cat "$prefix"*R2*.fastq.gz >"${prefix}_R2.fastq.gz"
done

或者，使用上面的sed代码，

printf '%s\n' *.fastq.gz | sed 's/^\([^_]*_[^_]*\).*/\1/' | uniq |
while read prefix; do
    cat "$prefix"*R1*.fastq.gz >"${prefix}_R1.fastq.gz"
    cat "$prefix"*R2*.fastq.gz >"${prefix}_R2.fastq.gz"
done

上面没有任何代码使用bash-specific (或特定于GNU的)特性，应该在所有POSIX中工作。

更新:我和生物信息工作者一起工作，我的一位同事评论道：

我们不应该简单地合并fastq文件..。在理想的世界中，人们应该分别映射每条车道，添加一个适当的RG，然后合并BAM文件。因为车道特定的影响是存在的，等等。它可以或多或少的重要，当然取决于下游的应用。

有关这方面的问题，请参考生物信息学栈交换站点。

Answer 2

Bash解决方案：

for f in *.fastq.gz; do 
    [[ "$f" =~ ^([^_]+_[^_]+)_.*(_[^_]+)_[0-9]+\.fastq\.gz$ ]]
    cat "$f" >> "${BASH_REMATCH[1]}${BASH_REMATCH[2]}.fastq.gz"
done

• ^([^_]+_[^_]+)_.*(_[^_]+)_[0-9]+\.fastq\.gz$ --将前2个前缀捕获到第一个捕获组(用于ex )的关键正则表达式。MC9_PREN和R-named后缀进入第二组被俘组(表示ex )。_R1)