Help

本节中我们提供一些和大作业的编程实现相关的补充材料。

0) 编程环境

• 由于完成本次作业需要一定的计算资源支持，我们为各个小组提供了T集群账户。关于T集群的使用，请大家参考 Run jobs in a cluster [Advanced]

• 完成大作业需要大家自己编写一些python，R和shell脚本。我们推荐大家根据目的的不同选择合适的编程环境。

• 有些脚本需要一定的计算量，并且会反复使用。这类脚本一经确定，一般就不会反复改动，而是会通过更换命令行参数改变输入输出以及脚本的行为。在实际计算中，运行这类脚本就不太适合使用交互式的工具，建议大家在进行实际计算前做好测试，在计算时用好python和R的命令行参数。

• 画图本身一般对计算资源要求较低，并且很可能需要反复更换不同的展示方法，尝试不同的参数等，所以在本地用jupyter notebook或rstudio这样的交互式的IDE会比较方便。建议大家把在服务器上计算好的，可以用来画图的结果以文本文档的形式保存，从而将有一定计算量的脚本和画图的脚本独立开来。

• 有一些IDE，vscode和pycharm，提供了远程连接的功能。如果有同学习惯使用IDE，又想利用远程的计算资源，可以尝试这些工具。

• 在T集群上，我们已经装好了一些常用工具，位于/data/2022-bioinfo-shared/softwares。大家可以将需要的工具添加到自己的环境变量中。

• 我们提供了预装了一些常用package的python和R的miniconda环境，可以通过如下方式进入:

eval "$(/data/2022-bioinfo-shared/softwares/miniconda3/bin/conda shell.bash hook)"
conda activate python-env # 进入python环境
conda deactivate #退出当前环境
conda activate r-env # 进入R环境

• 如果不进入conda环境，也可以通过指定相应的解释器使用这些package，例如将/data/2022-bioinfo-shared/softwares/miniconda3/envs/python-env/bin和/data/2022-bioinfo-shared/softwares/miniconda3/envs/r-env/bin加入$PATH环境变量。

1) 数据获取

我们需要用到的两个数据集分别在/data/2022-bioinfo-shared/data/quiz-I的exRNA-long和exRNA-short两个子目录下。在使用中大家不需要把这些数据都复制到自己的目录下，为了方便在自己的目录下访问，只需要建立一个符号链接即可:

cd  {working directory}
ln -s /data/2022-bioinfo-shared/data/quiz-I/exRNA-long .
ln -s /data/2022-bioinfo-shared/data/quiz-I/exRNA-short .

2) 构建count matrix

long RNA

对于长RNA，如前所述，我们提供的bam文件在基因组坐标上，所以大家用常规RNA-seq reads counting的方法即可，请参考2.1.expression-matrix一节对featureCounts的介绍。由于这373个样本中有一部分使用了同一个试剂盒的不同版本，所以有一部分样本是forward stranded，有一部分是reverse stranded，请大家参考metadata.txt的library一列。

# for reverse strand
export PATH=/data/2022-bioinfo-shared/softwares/subread/subread-2.0.3-Linux-x86_64/bin:$PATH
featureCounts --countReadPairs -O -M -s 2 -p -t exon -g gene_id -a genome/gff/gencode.v38.annotation.gff3 -o {output prefix} {bam}
# replace -s 2 with -s 1 for forward strand libraries

接下来就需要大家把不同样本的counts merge到一个count matrix中。大家可以自己实现，也可以使用我们提供的脚本(/data/2022-bioinfo-shared/data/quiz-I/scripts):

scripts/summarize-table.py --indir counts --formatter '{}' --sample-ids sample_ids.txt --row-field 0 --row-name gene_id --first-line 2 --value-field 6 --fillna --output count.matrix.txt

我们还要求大家统计不同RNA类型在long RNA reads中所占的比例。大家可以参考我们提供的脚本:

# 按mRNA,lncRNA,snoRNA,snRNA,srpRNA,tRNA,YRNA,intron,pseudogene的优先级把每个fragment assign到特定RNA类型
scripts/reads-assignment.py --input exRNA-long/bam/CRC-2418277.bam --strandness reverse --output CRC-2418277.txt --bed-dir exRNA-long/genome/bed

small RNA

对于小RNA，所有样本用的都是forward stranded的strand specific protocol。每一个样本对应9个bam文件，每个bam文件对应一种RNA类型。对于miRNA和piRNA，我们提供了一个脚本用来count assign到每一个转录本上的reads:

scripts/count-transcripts.py --bam {bam} --counts {counts} --stats {stats} -s forward

大家可能会注意到，在统计结果中invalid-strand一项在不同样本中都等于0，这是因为我们在用bowtie2 mapping的时候加了--norc参数，所以不会有reads被map到reverse strand上。

接下来同样需要大家把不同样本的counts merge到一个count matrix中，这样对于miRNA和piRNA都能产生一个count matrix。大家也可以参考我们提供的脚本来进行数据整理:

scripts/summarize-table.py --formatter "{}.txt" --indir output/counts --sample-ids sample_ids.txt --value-field 1 --row-field 0 --row-name miRNA_id --output miRNA.counts.txt  --fillna

我们可以把这两个矩阵合并到到一起用于后续分析。

3) 矩阵处理

过滤低表达基因

• 过滤低表达基因可以有不同的策略，我们给出一个使用edgeR的例子，鼓励大家多做其他尝试:

library(edgeR)
y <- edgeR::DGEList(counts=count.matrix)
keep <- edgeR::filterByExpr(y)
y <- y[keep, , keep.lib.sizes=FALSE]

表达量数据的normalization

• 用edgeR可以很方便的对数据进行归一化。请大家注意，edgeR::cpm这个函数虽然名字叫做CPM，它给出的归一化结果实际上取决于edgeR::calcNormFactors中的"method"。如果没有edgeR::calcNormFactors这一步，edgeR::cpm输出的才是我们一般说的CPM(count per million)。我们鼓励大家尝试不同的归一化方法，如edgeR提供的RLE, UQ等，以及其他一些基于内参基因的normalization等。

y <- edgeR::calcNormFactors(y,method="TMM")
cpm.matrix <- edgeR::cpm(y)

• RLE plot是一种可视化归一化结果的常见方法。EDASeq这个package直接提供了一个plotRLE的函数，感兴趣的同学也可以自己实现:

library(EDASeq)
plotRLE(cpm.matrix)

• 有了normalize后的数据我们就可以用PCA/MDS/tSNE等方法进行可视化，用metadata.txt中给出的信息标记样本点的颜色，从而分析那一个因素主导了数据的variations。具体实现请大家参考降维可视化一节。请大家注意，常用package的输入一般是行(row)为样本，列(column)为特征，对应表达矩阵的转置。

对于可视化和机器学习而言，我们一般会对CPM一类的表达量取对数作为输入。我们常常会遇到在一些样本中CPM为0的基因，这时直接取对数是没有定义的。对此，一个常见的解决办法是对所有的表达量都先加上一个"pseudocount"，再取对数，例如将log10(1+CPM)的数值当做是logCPM。如果提供log = TRUE参数，edgeR::cpm会对结果取对数，prior.count参数决定了pseudocount的数值大小。

去除批次效应

有很多R package都提供了去除批次效应的方法。以下提供了两个例子:

• 使用RUVSeq package中提供的RUVg函数根据control gene list进行batch correction

suppressPackageStartupMessages(library(RUVSeq))
# load tmm normalized matrix in log scale
log.tmm.matrix <- read.csv("small/miRNA.TMM.logged.txt",sep="\t",row.names = 1)
log.tmm.matrix <- as.matrix(log.tmm.matrix)
# calculate coefficient of variation
cv.index <- apply(log.tmm.matrix,1,function(x){sd(x)/mean(x)})
cv.index.sorted <- sort(cv.index)
# using 20% miRNA with smallest cv as empirical control gene as there is little knowledge for stably expressed miRNA
# for long RNA, your can use known housekeeping genes
empirical.control <- head(names(cv.index.sorted),as.integer(length(cv.index.sorted)*0.2))
# perform batch  correction with RUVg
res <- RUVg(log.tmm.matrix, empirical.control, k=2,isLog=TRUE)
log.tmm.matrix.ruvg <- res$normalizedCounts
log.tmm.matrix.ruvg <- as.data.frame(log.tmm.matrix.ruvg)
# saving results
write.table(log.tmm.matrix.ruvg,"small/miRNA.TMM.logged.ruvg.txt",sep="\t",quote = F)

• 使用sva package中提供的ComBat函数根据已知的batch信息进行batch correction

library(sva)
# load tmm normalized matrix in log scale
log.tmm.matrix <- read.csv("small/miRNA.TMM.logged.txt",sep="\t",row.names = 1)
log.tmm.matrix <- as.matrix(log.tmm.matrix)
metadata <- read.table("small/metadata.txt",sep="\t",row.names = 1,header = T)
log.tmm.matrix <- log.tmm.matrix[,row.names(metadata)]

# perform batch correction with ComBat
# only one batch allowed
log.tmm.matrix.combat <- ComBat(log.tmm.matrix, batch=metadata$library.prepration.day)
write.table(log.tmm.matrix.ruvg,"small/miRNA.TMM.logged.combat.txt",sep="\t",quote = F)

我们可以通过降维可视化分析降维可视化的效果。如果去除批次效应之前的样本点是按样本来源或实验批次分群的，经过这样的处理后，理想情况下，不同来源或批次的样本点会混在一起。下图给出了小RNA数据中miRNA表达量的例子。

对特征进行scaling

如果我们准备使用树模型之外的其他分类器，就需要先对特征进行scaling，再输入到模型中。请参考machine-learning-basics一节的相关介绍。

4) 特征选择

请大家参考feature engineering一节对特征选择的介绍。

请大家注意，特征选择原则上和模型拟合一样，不能用到测试集的label。

5) 模型评估和特征解释

请大家参考1.5.model-evaluation

6) peak calling

• 对于除了miRNA和piRNA之外的其他RNA类型，对于每个样本的每个bam文件，我们可以先用samtools去除duplication

#!/bin/bash
for sample_id in $(cat sample_ids.txt);do
  mkdir -p bam.nodup/$sample_id
  echo $sample_id
  for RNA in lncRNA mRNA snoRNA snRNA srpRNA tRNA YRNA;do
   samtools rmdup -s bam/$sample_id/${RNA}.bam bam.nodup/$sample_id/${RNA}.bam
  done
done

• 再将同一种RNA类型不同样本对应的bam文件merge到一起，并转换为bed文件:

#以lncRNA为例
ls  bam.nodup/*/lncRNA.bam > lncRNA.bamlist.txt
samtools merge -b lncRNA.bamlist.txt lncRNA.bam
bedtools bamtobed -i lncRNA.bam  > lncRNA.tags.bed

• 之后我们用piranha这个工具(安装路径:/data/2022-bioinfo-shared/softwares/piranha/piranha-1.2.1/bin)进行peak calling

• 将该工具的executable和它依赖的动态库加入环境变量:

export PATH=$PATH:/data/2022-bioinfo-shared/softwares/piranha/piranha-1.2.1/bin
export LD_LIBRARY_PATH=/data/2022-bioinfo-shared/softwares/gsl/lib:$LD_LIBRARY_PATH

• peak calling:

#以lncRNA为例
Piranha -b 20 -s lncRNA.tags.bed > lncRNA.peaks.bed

以上peak calling的过程供大家参考。另一种可行的方式是对每一个去了duplication的bam文件分别进行peak calling，再统计在多个样本中都可以检测到的peaks(即一些"recurrently identified peaks")，来定义一些"domain"。

• 我们把bam文件和计算出的peak load到IGV上进行可视化，进而对peak的形状有一些直观的认识

• 根据peak calling的结果，我们再用去duplication之前的bam文件统计落在每一个peak内部的reads数量。以下是一个样本lncRNA的例子:

bedtools coverage -a lncRNA.peaks.bed -b bam/Sample_1S1/lncRNA.bam -c > Sample_1S1.lncRNA.peaks.counts.txt

• 接下来对不同类型的RNA，同样需要大家将不同样本的counts合并到一个矩阵中，再将不同RNA类型对应的矩阵合并为一个矩阵。

6) Count reads for customized regions

• 定义read counting的区域有比较大的灵活性，大家可以直接下载已有的注释(注意应当在hg38的坐标上)或者自行定义一些区域

• 统计read数量的方法请大家参考我们之前对bedtools的介绍: 1.2 bedtools and samtools