Bowtie2比对过程中发生的一个小错误

测试开头

Bowtie2比对过程中发生的一个小错误今天是生信星球陪你的第629天

大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！

豆豆写于2020.5.21
记录一下踩的坑

上报错

手上有一些ATAC-seq的双端测序数据（25bp）进行了fastp质控过滤后准备进行比对，结果比对过程发生报错，信息如下

# Bowtie2
Error, fewer reads in file specified with -1 than in file specified with -2
terminate called after throwing an instance of 'int'
(ERR): bowtie2-align died with signal 6 (ABRT) (core dumped)

貌似显示双端测序数据中两个文件reads数不等

如何查看测序数据的信息

简单使用脚本

zcat  raw/test_1.fastq.gz | wc -l|  awk '{print $1/4}'
# 4237351
zcat  raw/test_2.fastq.gz | wc -l|  awk '{print $1/4}'
# 4237351

使用小工具

seqkit stat raw/test_1.fastq.gz raw/test_2.fastq.gz
# file               format  type   num_seqs      sum_len  min_len  avg_len  max_len
# raw/test_1.fastq.gz  FASTQ   DNA   4,237,351  105,933,775       25       25       25
# raw/test_2.fastq.gz  FASTQ   DNA   4,237,351  105,933,775       25       25       25

看到原始数据的两个文件中reads数相等

既然原始数据一致，那么质控过滤后的呢？

果然，利用seqkit stat发现clean的一对PE数据中reads数不一致

看了一下如何进行过滤的

sample=test
fq1=~/atac/raw/${sample}_1.fastq.gz
fastp -i $fq1 -o clean/${sample}.fq1.gz 
          -I $fq2 -O clean/${sample}.fq2.gz 
              --thread=4 --length_required=10 --n_base_limit=5 
              --compression=6 -R ${sample}_1  1>log/log.fastp.${sample}.1.txt 2>&1

fq2=~/atac/raw/${sample}_2.fastq.gz
fastp -i $fq2 -o clean/${sample}.fq2.gz 
              --thread=4 --length_required=10 --n_base_limit=5 
              --compression=6 -R ${sample}_2  1>log/log.fastp.${sample}.2.txt 2>&1

忽然想起来，这个脚本是从ChIP-seq脚本中复制粘贴过来的，只改了对应的fq1和fq2，但一般ChIP-seq都是单端的，所以这个fastp就是对单端数据的过滤

问题的根源在于：这里把双端数据当成了两个单端分别进行处理，那么结果就不一致啦

将代码修改为双端对应的就好了

# 如果是双端，就需要指定-i -o (for input 1) 
# -I -O (for input 2)
sample=test
fq1=~/atac/raw/${sample}_1.fastq.gz
fq2=~/atac/raw/${sample}_2.fastq.gz
fastp -i $fq1 -o clean/${sample}.fq1.gz 
      -I $fq2 -O clean/${sample}.fq2.gz 
              --thread=4 --length_required=10 --n_base_limit=5 
              --compression=6 -R ${sample}  1>log/log.fastp.${sample}.txt 2>&1

收获

一般来说，软件的报错信息都会显示在前面，需要观察加发散思维

例如

# bowtie2
Error: reads file does not look like a FASTQ file
terminate called after throwing an instance of 'int'
(ERR): bowtie2-align died with signal 6 (ABRT) (core dumped)

看到这里，就要去看原始的fq文件（主要信息就是第2行和第4行），可能存在长度不一致的情况

zcat fastq.gz | paste - - - - | awk -F"t" '{ if (length($2) != length($4)) print $0 }'
# 如果要挑出来这些reads
zcat fastq.gz | paste - - - - | awk -F"t" '{ if (length($2) != length($4)) print $0 }' |  tr 't' 'n' > error_reads.fastq