参照qczhang的github代码整理:github
在可以提交任务的节点(loginview02/mgt01)上,导入文件中的shell函数,然后像下面这样调用calc_icSHAPE这个函数。可以用bjobs来查看任务是否顺利提交。
OUT_DIR=$OUT_ROOT_Dir/mes_ch_vitro_yy
SAM_D1=/Share/home/zhangqf8/sunlei/data/icshape/mes/ch_vitro_transcript_low/Chr_D1.sam
SAM_D2=/Share/home/zhangqf8/sunlei/data/icshape/mes/ch_vitro_transcript_low/Chr_D2.sam
SAM_N1=/Share/home/zhangqf8/sunlei/data/icshape/mes/ch_vitro_transcript_low/Chr_vitro_N1.sam
SAM_N2=/Share/home/zhangqf8/sunlei/data/icshape/mes/ch_vitro_transcript_low/Chr_vitro_N2.sam
calc_icSHAPE $OUT_DIR $SAM_D1 $SAM_D2 $SAM_N1 $SAM_N2
readCollapse
参数:
-U $seqFile # 序列文件,用:分割
-o $seqCollapsed # 输出文件
-f $seqDatFasta # 把结果也存在一个自己定义的fasta文件中
去除所有读段中的重复读段。
首先,如果fastq文件太大,就用splitByBarcode把文件切成一个个部分;然后用readCollapse去除所有读段中重复出现的读段,这里所谓的重复读段就是指不管读段的质量怎么样,只要读段序列是一样的,就只保留一次。
readCollapse:输入文件是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。
输出的文件为四列:
计算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文件合并成一个文件,有下面的几种方法:
如果我们定义了头部和尾部的碱基要跳过多少个碱基的话,那些碱基就不参与上面的计算,直接设置为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。
对于每一条转录本,本来应该有大致12.5%的位点为1,但是发现很多转录本的最大值达不到1,这是因为那些RT值较大的位点多分布在转录本的前部分,然而这些位点的BD值却达不到200的阈值,导致这些位点被过滤为NULL。
这里的一个例子是把HEK293细胞胞质组分中的reads回帖到18S rRNA上
git clone https://github.com/lipan6461188/icSHAPE_protocol.git
cd icSHAPE_protocol
./icSHAPE_pipeline.sh run
# 把icSHAPE转成VARNA接受的格式
./icSHAPE_pipeline.sh cov
# 如果要清除空间
[./icSHAPE_pipeline.sh clear]
在图形界面下双击打开tools/VARNAv3-93.jar,把上面./icSHAPE_pipeline.sh cov
那一步的序列和结构复制粘贴进去,然后在窗口中右击选择下面的选项调整RNA的结构图形:
选择Color map => Load Values,输入文件rRNA.varna的路径,可以得到上色的18sRNA的结构。