目前已经是我接触生物信息学的第四年了，这里的内容由于学业压力也很久没有更新，现在可以趁着寒假的空闲时间写一写最近遇到的东西。我通常处理的都是RNA-seq数据，也全是来自illumina测序平台（自认为擅长RNA测序数据的个性化分析），但是最近学习ChIP-seq的时候，碰巧下载了一组SOLiD colorspace format格式的数据。由于一开始没有注意到这个问题，导致后面遇到了一系列的问题，非常多的报错，这里记录一下我的解决方案。
一、下载数据
这是一组果蝇的ChIP-seq数据，只下载了6个原始文件，包括野生型和突变型果蝇的两种抗体的ChIP，以及各自的Input，非常简单。使用aspera在ENA数据库下载数据。
$ ascp -l 100M -P 33001 -QT -k 1 -i
/这里是conda安装的位置/etc/asperaweb_id_dsa.openssh
era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR144/008/SRR1449908/SRR1449908.fastq.gz ./01_rawdata &
$ md5sum ./01_rawdata/*.fastq.gz
二、数据质量检测
使用fastqc和multiqc进行质量检测。
$ fastqc ./01_rawdata/*.fastq.gz -t 4 -o ./03_fastqc -f fastq &
$ multiqc ./03_fastqc -o ./03_multiqc &
到这一步都是一切正常，只是质量看上去...有点低，不过也能理解，毕竟是好多年前的数据了。
可以注意到序列中存在SOLID Adapter，至于这个adapter是什么，我也不知道。
三、数据过滤
本来想使用trim_galore进行过滤。
$ trim_galore -j 4 -q 20 --phred33 --length 25 -o ./04_cleandata
./01_rawdata/*.fastq.gz &
但是报错了。
File seems to be in SOLiD colorspace format which is not supported by Trim Galore (sequence is: 'T20103321021332323020223030003323033200023300033032')! Please use Cutadapt on colorspace files separately and check its documentation!
赶紧回去看原始的fastq文件，结果使人震惊：
啊这...
这种诡异的序列还是第一次遇到。在网上查了查，原来是colorspace数据，和我们通常看到的basespace不同。并且查到了一篇网上的文章：
（原来别人也遇到过，遇到这么少见的数据类型是不是说明我也算是身经百战了？）
现在我们换成cutadapt进行数据过滤。其中有一部分专门的colorspace选项：
-c, --colorspace：Enable colorspace mode: Also trim the color that is adjacent to the found adapter
--maq/--bwa：MAQ- and BWA-compatible colorspace output. This enables -c, -d, -t, --strip-f3 and -y '/1'. Enables colorspace mode (-c), double-encoding (-d), primer trimming (-t)
--no-zero-cap：不要将colorspace数据中的负质量值更改为零。colorspace读数的质量值有时是负的，BWA和Bowtie无法处理，默认情况下Cutadapt将colorspace数据中的负质量值转换为0。
--format=sra-fastq：当使用来自sra-toolkit包的FASTQ-dump程序将这些.sra文件转换为FASTQ格式时，colorspace序列将在每次开始时获得额外的质量值，使质量值数目比序列数多1。为了让cutadapt忽略额外的质量，在命令行中添加--format=sra-fastq。
我尝试了好几次cutadapt，最后确定了具体的参数。
$ cutadapt --bwa -a 330201030313112312 -q 20 -m 25 --max-n 2
--format=sra-fastq 01_rawdata/SRR1449908.fastq.gz
-o 04_cleandata/SRR1449908_clean.fastq.gz &
开始的时候我没有加-a，因为不知道adapter是什么。第一次尝试发现colorspace不支持-j。第二次用了-c和-t选项，而没有用--bwa。实际证明也是可以的。通过这一步之后再fastqc，确定了接头的序列是330201030313112312。然后重新cutadapt，并加上-a。但是由于我没有发现它的碱基质量行似乎多了一位，导致后面运行bowtie的时候报错，查看文档后加上了--format=sra-fastq。过滤后的数据再跑一次fastqc和multiqc。经过来来回回好几次cutadapt，才得到了看上去似乎可以的过滤后数据。
文档链接：https://cutadapt.readthedocs.io/en/v1.18/colorspace.html
四、比对
从网上查到bowtie支持colorspace数据，那我们就用bowtie进行比对吧。
# 建立索引
$ mkdir -p 02_index/bowtie_dm
$ bowtie-build --threads 4 -f
00_annotation/Drosophila_melanogaster.BDGP6.32.dna.toplevel.fa
02_index/bowtie_dm/bowtie_dm &
# 比对
$ gzip -d 04_cleandata/*.gz &
$ bowtie -p 4 -x 02_index/bowtie_dm/bowtie_dm
./04_cleandata/SRR1449909_clean.fastq -S ./06_bowtie/SRR1449909.sam
2>>./log/bowtie.log &
比对完一看，reads that failed to align 99.98%...
我又去网上查，发现bowtie处理colorspace数据需要加-C/--color，索引也是需要单独构建的。但是，对bowtie-build加上-C选项后，再次报错了。打印出了Usage，好像没有-C这个参数。这就有点使人疑惑了。我猜可能是版本问题，这个1.3.1版本可能不支持-C。我查看了bowtie的更新日志，找到了最后一次提到-C的版本是0.12.3。conda似乎无法给bowtie降级到0.12.3，所以我直接从浏览器把它下载下来了。
$ mkdir -p 02_index/bowtie_c
$ ./bowtie-0.12.3/bowtie-build -C -f
00_annotation/Drosophila_melanogaster.BDGP6.32.dna.toplevel.fa
02_index/bowtie_c/bowtie_c &
$ ./bowtie-0.12.3/bowtie -p 4 -C 02_index/bowtie_c/bowtie_c
./04_cleandata/SRR1449908_clean.fastq -S ./06_bowtie/SRR1449908.sam
2>>./log/bowtie.log &
最后得到了结果，reads with at least one reported alignment: 18051482 (72.09%)，好像也还行。接下来就是变成bam，然后拿去跑macs2了。
$ ls ./06_bowtie/*.sam | while read id ; do (samtools sort -@ 4
-O bam -o ./06_bowtie/$(basename $id ".sam").bam $id) ; done &
网友评论