从这一篇开始，是我重复一篇文献的记录，我在第五周文献学习笔记中提到过。
1. 过滤adapter
在第二篇笔记中，我判断接头序列是从raw data中看的，虽然准确但花了一些时间，而且如果待分析数据来自不同的建库，adapter序列就可能不一样了，这时候再看raw data确定接头，在时间上就不划算了。当时我就想，鉴于小RNA测序的特征（测序长度有35，36，49，50，51bp等，每条正常测序的reads总会包含接头），所以通过raw reads之间的比对，确定最保守的区域，不就是接头了吗？后来看到有推文就是这样做的。而MUSCLE比对工具使用，我之前就写过，所以这次确定接头就选择它了。
我选择了raw data中的50条序列（不要选最前面的几十条序列，测序质量不好），将其输入在线版的MUSCLE，输出格式为HTML，得到如下结果：
而具体trim的时候，adapter.fasta选前几个就够了，这次我选的：CTGTAGGCAC。仍然用的trimmomatic，但又有了新的认识，关于其单端过滤模式的ILLUMINACLIP选项。
#1
#ILLUMINACLIP:SE.fa:1:30:6
java -jar trimmomatic-0.36.jar SE -threads 4 -phred33 SRR3269248.fastq.gz out.SRR3269248.fastq.gz ILLUMINACLIP:SE.fa:1:30:6
#2
#ILLUMINACLIP:SE.fa:2:30:10
java -jar trimmomatic-0.36.jar SE -threads 4 -phred33 SRR3269248.fastq.gz out2.SRR3269248.fastq.gz ILLUMINACLIP:SE.fa:2:30:10
#3 最终选择
#ILLUMINACLIP:SE.fa:1:30:6 SLIDINGWINDOW:5:20 MINLEN:18 # 5:20 表示在5nt的滑动窗口中平均质量20以下的窗口会被切去
java -jar trimmomatic-0.36.jar SE -threads 4 -phred33 SRR3269248.fastq.gz out3.SRR3269248.fastq.gz ILLUMINACLIP:SE.fa:1:30:6 SLIDINGWINDOW:5:20 MINLEN:18
以第1次尝试为例，ILLUMINACLIP的选项中的1表示允许adapter与reads比对时有1个错配，30这个值（不一定就是30）只在双端测序中起作用，但单端模式也要加上，6表示最小匹配分值，它是这样计算的：raw reads的长度为36，结合小RNA的长度范围，可知一个raw reads带有adapter的长度在10左右到18之间，匹配一个碱基得0.6分，(10-1)*0.6，所以这个值我定的是6。这也是为什么1，3图的长度分布符合预期，在21，24处出现峰值。反观第二次尝试，第三个值定的10，基本不可能符合，所以最后发现adapter仍然在clean.reads文件中大量存在，而我在网上找的很多教程包括官方手册都没有指出这个SE模式下的陷阱。1，3图的区别不大，加了测序质量控制，reads末端若质量低会被切去。
2. reads长度分布图
我在一些文献中看到过这种做法，其目的就是比较不同组织不同状态下的各样本之间所表达的小RNA的长度变化。这篇文献的补充材料图2就是这个图，还是3维的。
#我的命令行如下
nohup ls out.SRR3269*.fastq.gz \
| while read id; do less ${id} \
| awk '{ if(NR%4==2) print length($0) }' \
| sort -n | uniq -c \
| awk '{ OFS="\t";$1=$1;print $2"\t"$1 }' > len.${id%.*}.txt; done &
在这之后就能得到每个样本的reads长度分布表格了，
以为可以开始画图了，结果还是太年轻。我居然看到了在27nt出现峰值的文件（raw reads长度为42），理论上不是21或24嘛。我赶紧又找了一下这些42读长的fastq的接头。吓了一跳，原来还有barcode。
样本与样本之间不一样，所以其作用应该是在一次上机测序中区分样本的。故重新去barcode之后去接头。
#加上 HEADCROP:3
java -jar trimmomatic-0.36.jar SE -threads 4 -phred33 xxx.fastq.gz out.xxx.fastq.gz ILLUMINACLIP:SE.fa:1:30:6 SLIDINGWINDOW:5:20 MINLEN:18 HEADCROP:3
继续整理。
#先求每个样本clean_reads总数
ls len.out.SRR3269* | while read id; do awk '{sum+=$2}END{print sum}' $id; done | less
#每个样本18-28nt的reads数
$ ls len.out.SRR3269* | while read id
> do
> less $id | awk '{ if ($1>=18 && $1<=28) print $0 }' > 18_28.${id}
> done
#合并
paste 18_28.len.out.SRR3269* > tmp_59sample_18_28_reads_count.txt
#提取数据列
awk '{ for(i=2;i<=118;i=i+2) printf("%s\t",$i);printf("\n") }' tmp_59sample_18_28_reads_count.txt > 59sample_18_28_reads_count.txt
简单整理之后得到如下表格，接着导入R中做进一步处理，折线图就行了，这个3D图我暂时还不会画。
3. bowtie比对reads到参考基因组OR已知miRNA/MIRNA序列
Hairpin前体序列和miRNA成熟序列都可以在miRBase数据库下载gff后提取得到。因为gff中包含+-链的信息，所以我用两个模式提取，后面再来比较区别。
#可能需要先改个名字
#sed -i 's/>/>chr/g' Arabidopsis_thaliana.TAIR10.dna.toplevel.fa
#不加-s全部都取参考基因组上面的序列（+链）
bedtools getfasta -fi Arabidopsis_thaliana.TAIR10.dna.toplevel.fa -bed ath_MIRNA.gff3 -fo ath_MIRNA.fasta
eg:
less ath_MIRNA.fasta | grep ">" | head -n 2
>chr1:28499-28706
>chr1:78926-79037
#加了-s之后，对于-链的gff记录会取反向互补序列
bedtools getfasta -s -fi Arabidopsis_thaliana.TAIR10.dna.toplevel.fa -bed ath_MIRNA.gff3 -fo ath_MIRNA+-.fasta
eg:
less ath_MIRNA+-.fasta | grep ">" | head -n 2
>chr1:28499-28706(+)
>chr1:78926-79037(-)
接下来将小RNA测序reads比对到miRBase数据库上面的前体序列。
#建立索引
bowtie-build Arabidopsis_thaliana.TAIR10.dna.toplevel.fa Arabidopsis_thaliana.TAIR10.dna.toplevel
bowtie-build ath_MIRNA.fasta ath_MIRNA
bowtie-build ath_MIRNA+-.fasta ath_MIRNA+-
#比对到前体序列
1. bowtie -v 1 -p 4 -m 6 --norc ath_MIRNA out.SRR3269248.fastq.gz SRR3269248.bwt
2. bowtie -v 1 -p 4 -m 6 --norc ath_MIRNA+- out.SRR3269248.fastq.gz SRR3269248.2.bwt
参数解释：
-v: 允许一个错配
-p: 四个线程
-m: 当多比对超过这个数时，认为是未必对
--norc: 仅对提供的参考序列做比对，不会对参考序列求负链后再跟负链比对，此处需加上，这样上面两种比对才有区别。
前面下载前体序列后数了一下，共326条前体序列。现在对上面的比对结果简单统计一下。
#多少个已知的前体序列在此次比对中被覆盖了
less SRR3269248.bwt | cut -f3 | sort | uniq -c | wc -l
228
less SRR3269248.2.bwt | cut -f3 | sort | uniq -c | wc -l
283
可以看出，按照gff中提供的+-链信息提取序列后做参考，比简单地全部按照+链提取准确全面。
前面提到了--norc这个参数，那如果我们允许参考序列的+-链都参与比对呢，结果会怎么样？
bowtie -v 1 -p 4 -m 6 ath_MIRNA out.SRR3269248.fastq.gz SRR3269248.3.bwt
less SRR3269248.3.bwt | cut -f3 | sort | uniq -c | wc -l
284
和上面的第二种比对接近。如果我们只知道产生miRNA的区域，并不知道它是参考基因组的+/-链，也许定量的时候直接提取+链，再正负链比对，应该也可以很好地得出定量结果。
说到定量当然要求出落在每一个参考前体序列上的reads数，先看看bowtie的输出文件。
下面是第三种比对得到的结果，只取了两个记录。比对到正链的结果较好识别，着重看看比对到负链的
SRR5096869-5404581-1 + zma-miR2118e/.22 0 TTCCTGATGTCTCCCATTCCTA ?@@?;DDDH=CDFHGIGGIIGI 0
SRR5096869-17330461-1 - zma-miR2118e/.22 0 TTCCTGATGCCTCCCATT IIIJJHHHHHFFFFFBCC 0 8:T>C
#这里的8是怎么来的？
已知参考序列是成熟miRNA序列，ref.fasta文件中记录为： TTCCTGATG-T-CTCCCATTCCTA
bwt比对结果中的第5列： TTCCTGATG-C-CTCCCATT（5' -> 3'排列，第10个碱基变化，T -> C；并且是原始序列的反向互补序列；反向数第9位，也就编号为8）
原始小RNA的reads文件，fastq文件： AATGGGAG-G-CATCAGGAA（回到这条原始序列中，它仍然是5' -> 3'排列，第9个碱基变化，编号为8）
这两种理解都能确定是8
这两行命令都能求reads count，之后就能将表达量标准化了。
1. less SRR3269248.3.bwt | perl -alne '{$h{$F[2]}++}END{print "$_\t$h{$_}" foreach sort keys %h}' | cat
2. less SRR3269248.3.bwt | cut -f3 | sort | uniq -c | cat
以上只是求了一个样本，将每一个样本都这样做，就能知道哪些组织特异表达哪些miRNA了。然后将所有的样本中验证的miRBase中的miRNA取并集就能知道这个研究验证了miRBase中的多少个miRNA。
另外，上面只是对已知的miRNA定量和验证，不能发现新的miRNA。要想发现新的miRNA，需要对基因组比对。之前有一篇笔记简单地写过用软件套装来预测miRNA，见笔记4。
reference
bowtie输出格式详见：http://blog.sina.com.cn/s/blog_6d4b5abd0101if53.html
