目前已經(jīng)是我接觸生物信息學(xué)的第四年了稻轨，這里的內(nèi)容由于學(xué)業(yè)壓力也很久沒有更新线欲，現(xiàn)在可以趁著寒假的空閑時(shí)間寫一寫最近遇到的東西圾旨。我通常處理的都是RNA-seq數(shù)據(jù)，也全是來自illumina測序平臺(tái)（自認(rèn)為擅長RNA測序數(shù)據(jù)的個(gè)性化分析）谣旁，但是最近學(xué)習(xí)ChIP-seq的時(shí)候昔头，碰巧下載了一組SOLiD colorspace format格式的數(shù)據(jù)驮吱。由于一開始沒有注意到這個(gè)問題拇砰，導(dǎo)致后面遇到了一系列的問題，非常多的報(bào)錯(cuò)甥材，這里記錄一下我的解決方案芝发。
一、下載數(shù)據(jù)
這是一組果蠅的ChIP-seq數(shù)據(jù)悟耘，只下載了6個(gè)原始文件管行，包括野生型和突變型果蠅的兩種抗體的ChIP，以及各自的Input，非常簡單极颓。使用aspera在ENA數(shù)據(jù)庫下載數(shù)據(jù)闰围。
$ 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
二忠售、數(shù)據(jù)質(zhì)量檢測
使用fastqc和multiqc進(jìn)行質(zhì)量檢測。
$ fastqc ./01_rawdata/*.fastq.gz -t 4 -o ./03_fastqc -f fastq &
$ multiqc ./03_fastqc -o ./03_multiqc &
到這一步都是一切正常睬捶，只是質(zhì)量看上去...有點(diǎn)低积仗，不過也能理解，畢竟是好多年前的數(shù)據(jù)了。
可以注意到序列中存在SOLID Adapter贩耐，至于這個(gè)adapter是什么更鲁，我也不知道樊拓。
三赴叹、數(shù)據(jù)過濾
本來想使用trim_galore進(jìn)行過濾蚕冬。
$ trim_galore -j 4 -q 20 --phred33 --length 25 -o ./04_cleandata
./01_rawdata/*.fastq.gz &
但是報(bào)錯(cuò)了锨苏。
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文件驯击，結(jié)果使人震驚：
這種詭異的序列還是第一次遇到择吊。在網(wǎng)上查了查所森，原來是colorspace數(shù)據(jù)掩幢，和我們通常看到的basespace不同度硝。并且查到了一篇網(wǎng)上的文章：
現(xiàn)在我們換成cutadapt進(jìn)行數(shù)據(jù)過濾折晦。其中有一部分專門的colorspace選項(xiàng)：
-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數(shù)據(jù)中的負(fù)質(zhì)量值更改為零。colorspace讀數(shù)的質(zhì)量值有時(shí)是負(fù)的，BWA和Bowtie無法處理锅论，默認(rèn)情況下Cutadapt將colorspace數(shù)據(jù)中的負(fù)質(zhì)量值轉(zhuǎn)換為0剔猿。
--format=sra-fastq：當(dāng)使用來自sra-toolkit包的FASTQ-dump程序?qū)⑦@些.sra文件轉(zhuǎn)換為FASTQ格式時(shí)，colorspace序列將在每次開始時(shí)獲得額外的質(zhì)量值慌闭，使質(zhì)量值數(shù)目比序列數(shù)多1豺妓。為了讓cutadapt忽略額外的質(zhì)量，在命令行中添加--format=sra-fastq耕魄。
我嘗試了好幾次cutadapt蔼两，最后確定了具體的參數(shù)。
$ 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 &
開始的時(shí)候我沒有加-a渐北，因?yàn)椴恢繿dapter是什么破托。第一次嘗試發(fā)現(xiàn)colorspace不支持-j。第二次用了-c和-t選項(xiàng)敌呈，而沒有用--bwa。實(shí)際證明也是可以的。通過這一步之后再fastqc，確定了接頭的序列是330201030313112312卑吭。然后重新cutadapt唉窃，并加上-a件已。但是由于我沒有發(fā)現(xiàn)它的堿基質(zhì)量行似乎多了一位，導(dǎo)致后面運(yùn)行bowtie的時(shí)候報(bào)錯(cuò)连茧，查看文檔后加上了--format=sra-fastq虐唠。過濾后的數(shù)據(jù)再跑一次fastqc和multiqc溉愁。經(jīng)過來來回回好幾次cutadapt家肯，才得到了看上去似乎可以的過濾后數(shù)據(jù)娘汞。
文檔鏈接：https://cutadapt.readthedocs.io/en/v1.18/colorspace.html
四、比對
從網(wǎng)上查到bowtie支持colorspace數(shù)據(jù)狸捅，那我們就用bowtie進(jìn)行比對吧。
# 建立索引
$ 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%...
我又去網(wǎng)上查输莺，發(fā)現(xiàn)bowtie處理colorspace數(shù)據(jù)需要加-C/--color啰脚，索引也是需要單獨(dú)構(gòu)建的。但是露氮，對bowtie-build加上-C選項(xiàng)后三妈，再次報(bào)錯(cuò)了涧窒。打印出了Usage，好像沒有-C這個(gè)參數(shù)。這就有點(diǎn)使人疑惑了金闽。我猜可能是版本問題，這個(gè)1.3.1版本可能不支持-C淘菩。我查看了bowtie的更新日志，找到了最后一次提到-C的版本是0.12.3阿蝶。conda似乎無法給bowtie降級(jí)到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 &
最后得到了結(jié)果匕积，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 &
