icSHAPE数据处理流程

参照qczhang的github代码整理:github

在集群上并行地运行icSHAPE的流程

在可以提交任务的节点(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

输入数据

  1. 原始测序读段数据(DMSO的两个重复 + NAI的两个重复):Fastq
  2. 转录组:Fasta
  3. (可选)注释文件:GTF

处理流程

  1. readCollapse

    参数:

    -U $seqFile            # 序列文件,用:分割
    -o $seqCollapsed       # 输出文件
    -f $seqDatFasta        # 把结果也存在一个自己定义的fasta文件中

    ​ 去除所有读段中的重复读段。

    ​ 首先,如果fastq文件太大,就用splitByBarcode把文件切成一个个部分;然后用readCollapse去除所有读段中重复出现的读段,这里所谓的重复读段就是指不管读段的质量怎么样,只要读段序列是一样的,就只保留一次。

    ​ readCollapse:输入文件是fastq文件,输出文件是下面的格式:

  2. 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

           去除读段中的部分碱基,这主要分为两个步骤:

    当定义了一个文件中包含一些Adaptor读段时,就可以使用Trimmomatic把读段中得Adaptor去掉。

  3. 比对读段到转录组上

    Bowtie2 \
    -U $seqTrimmed \          # Unpaired Reads(fastq)
    -S $outputSamFile \       # 输出SAM文件
    -x ref_genome_index \     # 参考基因组索引
    --non-deterministic \     # 以当前的时间生成随机数种子。由于前面去除相同的reads,所以这个参数没有必要的
    --time              \     # 比对时输出时间信息
    --norc                    # 比对找转录组上时,强烈推荐使用这个参数

    得到的是一个SAM文件。当遇到一个Read可以同时map到多个位置时,Bowtie2在最好的Map中随机地挑选一个。

  4. 计算RPKM

    关于SAM文件的格式,可以参考:http://genome.sph.umich.edu/wiki/SAM

  5. 计算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的步骤中,会去除转录本第一位的信息,但是这一位的重要性就在于影响计算尺度因子。

  6. combineRTreplicates.pl

    参数为:

    -i $allInputSignalFile               # 要结合的文件,用:号分割
    -o $combinedInputSignalFile          # 输出文件

    把两个DMSO和两个NAI-N3的文件合并在一起。合并的方法就是简单的把对应的转录本上的RT和BD加在一起。如果某些转录本只在两次重复中的一次中出现,则也输出结果。如下图所示:

    上面的ENST00000395829只在其中的一个文件中出现了,则只有一个rpkm。从下面可以看出重复实验之间仅仅是简单地相加:

  7. 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)。计算尺度因子的方法:

    1. 把有效序列(>40bp)的值从大到小排序;
    2. 有下面的四种方法来得到一段采样序列:
      • smart:从大到小收集所有的非零数字;
      • upper:从大到小收集所有非零数字中的前半部分;
      • decile(n):取从大到小的数字中(n-1)/10 ~ n/10;
      • vigintile(n):取从大到小的数字中(n-1)/20 ~ n/20;

    比如上面的vigitile2就是取从大到小排序的非零数字中的1/20 ~ 2/20。

    1. 通过采样序列得到尺度因子,有下面的三种方法:
      • median:取采样序列中的中值;
      • mean:取采样序列中的均值;
      • peak:取采样序列中的最大值;

    ​ 得到尺度因子ScalFac后,把它除以scaling form:

    \(\$scalingFactor=\frac{\$scalingFactor}{\$scalingForm}\)

    下面把每一个值用ScallingFactor来校正:

    1. 如果normalizeMethod是log: \(x_i = log_2(\frac{x_i}{scalingFactor}+1)\)
    2. 没有定义log时 : \(x_i = \frac{x_i}{scalingFactor}\)
  8. 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\)

    下面做如下处理:

    1. 首先把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。

    2. 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
  9. 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值,步骤如下:

    1. 计算上面文件中第五列的第二项的所有非NULL值的均值,把这个值作为coverage:

    my $coverage = eval ( join ( "+",@fgrtstops ) ) / scalar ( @fgrtstops );

    1. 果coverage大于上面定义的-T2,则继续下面的步骤;
    2. 把上面第三列的几个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的结构。