readCollapse
参数：
-U $seqFile # 序列文件,用:分割
-o $seqCollapsed # 输出文件
-f $seqDatFasta # 把结果也存在一个自己定义的fasta文件中 去除所有读段中的重复读段。
首先，如果fastq文件太大，就用splitByBarcode把文件切成一个个部分；然后用readCollapse去除所有读段中重复出现的读段，这里所谓的重复读段就是指不管读段的质量怎么样，只要读段序列是一样的，就只保留一次。
readCollapse：输入文件是fastq文件，输出文件是下面的格式：
- 自己定义的fasta文件（两列）：第一列：读段出现的次数；第二列：序列
- fastq文件：保留fastq文件中第一次出现的fastq读段
Trimming
参数：
-U $seqCollapsed # 序列文件
-o $seqTrimmed # 输出trimmed序列文件
-l 13 # crop of leading nucleotides
-t 0 # crop of tailing nucleotides
-c phred33 # fastq coding
-a $ICSHAPE/data/TruSeq2-PE.fa # Adaptor文件位置
-m 25 # min length去除读段中的部分碱基，这主要分为两个步骤：
simpleTrim简单地去除两端的碱基
需要定义一个读段最小长度，被trim过的读段需要大于这个长度才能被保留。读段头部的部分碱基被武断地切除。
trimmomatic除去Adaptor
当定义了一个文件中包含一些Adaptor读段时，就可以使用Trimmomatic把读段中得Adaptor去掉。
比对读段到转录组上
Bowtie2 \
-U $seqTrimmed \ # Unpaired Reads(fastq)
-S $outputSamFile \ # 输出SAM文件
-x ref_genome_index \ # 参考基因组索引
--non-deterministic \ # 以当前的时间生成随机数种子。由于前面去除相同的reads，所以这个参数没有必要的
--time \ # 比对时输出时间信息
--norc # 比对找转录组上时，强烈推荐使用这个参数得到的是一个SAM文件。当遇到一个Read可以同时map到多个位置时，Bowtie2在最好的Map中随机地挑选一个。
计算RPKM
关于SAM文件的格式，可以参考：http://genome.sph.umich.edu/wiki/SAM
首先从@SQ开头的行中读取转录本的名字和长度，存在一个Hash中；
读取Read的Map信息，如果Tag是0或者16，表示这个Single-End测序的读段Map上去的，如果是99或355，就是Paired-End Map上去的，这里只考虑Single-End的情况（关于Tag的信息参考）
Tag为4表示Read没有Map上，所以不考虑，只考虑Tag为0或16的读段。如果读段是Single Map的，那只会在这个文件中出现一次；如果读段是Multiple Map的，那在文件中会同时连续出现（这一点特征很重要）。我们把读段的名字、Map上的转录本的名字都记录到转录本的Hash中，统计每一个转录本被Map上的读段数目。如果读段是Single Map的，则给转录本加1，如果读段同时map到3个不同的转录本上，则给每一个转录本加1/3。
计算每一个转录本的RPKM：
其中，100亿表示转录本的长度应该除以100万，总的Reads数应该除以1000。
输出的文件为四列：
- transcriptid； 2. length；3. uniquereadnum；4. multipleread_num；5. RPKM
计算RT值
参数为：
-i $samFile # Bowtie产生的SAM文件
-o $rtFile # 输出文件
-r $rpkmFile # 上一步产生的RPKM文件
-c $config{MINLOAD}=5 # RPKM的最小值 这一步主要是把每一条转录本上每一个碱基被Map到的次数作为该碱基的Base Density；而每一条Reads停在的位置作为RT。
比如上面的文件中，第一行是Base Density，有两条Reads从第二个碱基起被Map，所以可以推断出第一个碱基可能被NAI-N3修饰了两次，所以第一个碱基的RT是2。只有RPKM大于上面的-c参数5时才保留这条转录本的信息。
这里有一个值得注意的地方就是，生成的文件中后面的RT和BaseDensity的数量是转录本的长度+1，这一位就加在第一位的位置（是一个没有意义的碱基），这样第一位碱基的BD总是0，最后一位的碱基的RT总是0。在后面normalize的步骤中，会去除转录本第一位的信息，但是这一位的重要性就在于影响计算尺度因子。
combineRTreplicates.pl
参数为：
-i $allInputSignalFile # 要结合的文件，用:号分割
-o $combinedInputSignalFile # 输出文件把两个DMSO和两个NAI-N3的文件合并在一起。合并的方法就是简单的把对应的转录本上的RT和BD加在一起。如果某些转录本只在两次重复中的一次中出现，则也输出结果。如下图所示：
上面的ENST00000395829只在其中的一个文件中出现了，则只有一个rpkm。从下面可以看出重复实验之间仅仅是简单地相加：
Normalize RT
参数为：
-i $combinedInputSignalFile # DMSO File
-o $normalizedInputSignalFile # NAI File
-m mean:vigintile2 # normalize method
-d 32 # head to skip
-l 32 # tail to skip
-f 100 # scaling form (the benchmarked value will be scaled to this)(没用)输出的结果为：
要求除去Head部分和Tail部分的碱基后有效转录本的长度应该控制在40bp以上，如果不符合这个条件，则上面的-d和-l控制的参数就就会减半，直到满足条件为止。取40bp以上的有效长度的作用是用来计算尺度因子（ScallingFactor）。计算尺度因子的方法：
- 把有效序列（>40bp）的值从大到小排序；
- 有下面的四种方法来得到一段采样序列：
- smart：从大到小收集所有的非零数字；
- upper：从大到小收集所有非零数字中的前半部分；
- decile(n)：取从大到小的数字中(n-1)/10 ~ n/10；
- vigintile(n)：取从大到小的数字中(n-1)/20 ~ n/20；
比如上面的vigitile2就是取从大到小排序的非零数字中的1/20 ~ 2/20。
- 通过采样序列得到尺度因子，有下面的三种方法：
- median：取采样序列中的中值；
- mean：取采样序列中的均值；
- peak：取采样序列中的最大值；
得到尺度因子ScalFac后，把它除以scaling form：
下面把每一个值用ScallingFactor来校正:
- 如果normalizeMethod是log: \(x_i = log_2(\frac{x_i}{scalingFactor}+1)\)
- 没有定义log时 : \(x_i = \frac{x_i}{scalingFactor}\)
calcEnrich
参数：
-f $normalizedTargetSignalFile # NAI-N3文件
-b $normalizedInputSignalFile # DMSO文件
-o $icShapeAllFile # 输出文件
-w factor5:scaling1 # winsor method, defined as, e.g. factor5:scaling10
-y 10 # dividing factor
-x 0.25 # subtraction factor
-e complex # enrich method, subtraction or dividing
[-g 1] # 上一步是log归一化方法把其中的参数用字母表示：
wisor factor(5): U(95%), L(5%)
dividingfactor(10): \(\alpha_1\)
substractionfactor(0.25): \(\alpha_0\)
winsor scalling(1): \(\beta\)
下面做如下处理：
首先把DMSO文件和NAI-N3文件合并成一个文件，有下面的几种方法：
- 如果-g参数是1或者-e参数是subtraction，则把每一个位点做如下处理：
\(fg\_rt-\alpha_0*bg\_rt\) - 如果-e参数是dividing，则把每一个位点做如下处理：
\(\alpha_1*\frac{fg\_bd}{bg\_bd}\) - 如果-e参数是其他，则把每一个位点做如下处理：
\(\frac{\alpha_1*(fg\_rt-\alpha_0*bg\_rt)}{bg\_bd}\) 如果我们定义了头部和尾部的碱基要跳过多少个碱基的话，那些碱基就不参与上面的计算，直接设置为NULL。在上面的后面两点中，当bg_db为0时，直接把值设置为NULL。
winsorization
这一部分：
my $headToSkip = 5; my $tailToSkip = 32;winsorization首先需要计算出winsorwindow：
my ( $winsorUpper, $winsorLower ) = &winsorWindow ( $ref_transcript_signal->{$trans}, $headToSkip, $trimed_last, $winsorFactor=5)这里主要是这么做的：首先取去掉头部去掉尾部的末端序列，然后把所有的非NULL值从大到小排序，取排在其中的第winsorFactor/100的值作为Upper(U)，排在其中的第1-winsorFactor/100的值作为Lower(L)
控制L的值大于等于0，所有的值中大于U的全部设置为U，小于L的设置为L，再对所有的非NULL碱基如果定义了scalling(1)，就做如下处理：
\(\beta\frac{x_i-L}{U-L}\)
头部和尾部跳过的碱基不参与计算winsorWindow，但是所有的非NULL碱基都应该进行上面的操作。
输出的结果如下所示：
第一列 trans_id 第二列 trans_length 第三列 fg_rpkm1 fg_rpkm2 bg_rpkm1 bg_rpkm2 第四列 fg_rt_ScalFact( A ) bg_bd_ScalFact( B ) bg_rt_ScalFact( C ) 第五列 base_enrichment A*fg_rt C*bg_rt B*bg_bd
FilterEnrichment
参数：
-i $icShapeAllFile # 上一步的输入文件
-o $icShapeFile # 输出文件
-t 200 # threshold to select window for correlation calculation
-T 2 # threshold of coverage
-s 5 # head to skip
-e 30 # end to skip这一步主要是过滤掉一些bg_BD过小的位点的RT值，步骤如下：
- 计算上面文件中第五列的第二项的所有非NULL值的均值，把这个值作为coverage:
my $coverage = eval ( join ( "+",@fgrtstops ) ) / scalar ( @fgrtstops );
- 果coverage大于上面定义的-T2，则继续下面的步骤；
- 把上面第三列的几个RPKM取均值，作为最终结果中的rpkm；对于每一个非头非尾的（5<\(X_i\)<len-30）碱基，如果它的第五列中B*bg_bd那一项大于-t设置的阈值，则输出，否则输出NULL。
