# Help

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

## 0) 编程环境

* 由于完成本次作业需要一定的计算资源支持，我们为各个小组提供了T集群账户。关于T集群的使用，请大家参考 [cluster](https://book.ncrnalab.org/teaching/getting-started/cluster "mention")
* 完成大作业需要大家自己编写一些python，R和shell脚本。我们推荐大家根据目的的不同选择合适的编程环境。
* 有些脚本需要一定的计算量，并且会反复使用。这类脚本一经确定，一般就不会反复改动，而是会通过更换命令行参数改变输入输出以及脚本的行为。在实际计算中，运行这类脚本就不太适合使用交互式的工具，建议大家在进行实际计算前做好测试，在计算时用好python和R的命令行参数。
* 画图本身一般对计算资源要求较低，并且很可能需要反复更换不同的展示方法，尝试不同的参数等，所以在本地用jupyter notebook或rstudio这样的交互式的IDE会比较方便。建议大家把在服务器上计算好的，可以用来画图的结果以文本文档的形式保存，从而将有一定计算量的脚本和画图的脚本独立开来。
* 有一些IDE，vscode和pycharm，提供了远程连接的功能。如果有同学习惯使用IDE，又想利用远程的计算资源，可以尝试这些工具。
* 在T集群上，我们已经装好了一些常用工具，位于`/data/2022-bioinfo-shared/softwares`。大家可以将需要的工具添加到自己的环境变量中。
* 我们提供了预装了一些常用package的python和R的miniconda环境，可以通过如下方式进入:

```bash
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`两个子目录下。在使用中大家**不需要**把这些数据都复制到自己的目录下，为了方便在自己的目录下访问，只需要建立一个符号链接即可:

```bash
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](https://book.ncrnalab.org/teaching/part-iii.-ngs-data-analyses/2.rna-seq/2.1.expression-matrix)一节对featureCounts的介绍。由于这373个样本中有一部分使用了同一个试剂盒的不同版本，所以有一部分样本是forward stranded，有一部分是reverse stranded，请大家参考`metadata.txt`的`library`一列。

```bash
# 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`):

```bash
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中所占的比例。大家可以参考我们提供的脚本:

```bash
# 按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:

```bash
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。大家也可以参考我们提供的脚本来进行数据整理:

```bash
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的例子，鼓励大家多做其他尝试:

```r
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等。

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

* RLE plot是一种可视化归一化结果的常见方法。[EDASeq](https://bioconductor.org/packages/EDASeq)这个package直接提供了一个plotRLE的函数，感兴趣的同学也可以自己实现:

```r
library(EDASeq)
plotRLE(cpm.matrix)
```

![](https://4115668567-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-LPVsf5VZbQ7h14X29qW%2Fuploads%2Fgit-blob-5a12efc717151ed3e77ae874e40f3196c5f3940f%2Fquiz.RLE-plot.png?alt=media)

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

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

### 去除批次效应

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

* 使用[RUVSeq](https://bioconductor.org/packages/RUVSeq) package中提供的RUVg函数根据control gene list进行batch correction

```r
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](https://bioconductor.org/packages/release/bioc/html/sva.html) package中提供的ComBat函数根据已知的batch信息进行batch correction

```r
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表达量的例子。

![](https://4115668567-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-LPVsf5VZbQ7h14X29qW%2Fuploads%2Fgit-blob-61a8d4748c2893b240152fd36a7408d6f237ec88%2Fquiz.batch.correction.PCA.png?alt=media)

### 对特征进行scaling

如果我们准备使用树模型之外的其他分类器，就需要先对特征进行scaling，再输入到模型中。请参考[machine-learning-basics](https://book.ncrnalab.org/teaching/part-iv.-machine-learning/1.machine-learning-basics/1.1.data-preprocessing#2-shu-ju-de-scaling)一节的相关介绍。

## 4) 特征选择

请大家参考[feature engineering](https://book.ncrnalab.org/teaching/part-iv.-machine-learning/1.machine-learning-basics/1.3.feature-engineering)一节对特征选择的介绍。

{% hint style="warning" %}
请大家注意，特征选择原则上和模型拟合一样，不能用到测试集的label。
{% endhint %}

## 5) 模型评估和特征解释

请大家参考[1.5.model-evaluation](https://book.ncrnalab.org/teaching/part-iv.-machine-learning/1.machine-learning-basics/1.5.model-evaluation)

## 6) peak calling

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

```bash
#!/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文件:

```bash
#以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:

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

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

* 我们把bam文件和计算出的peak load到IGV上进行可视化，进而对peak的形状有一些直观的认识
* 根据peak calling的结果，我们再用去duplication之前的bam文件统计落在每一个peak内部的reads数量。以下是一个样本lncRNA的例子:

```bash
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-samtools](https://book.ncrnalab.org/teaching/part-iii.-ngs-data-analyses/1.mapping/1.2-bedtools-samtools "mention")