# Help

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

## 0) 编程环境

* 由于完成本次作业需要一定的计算资源支持，我们为各个小组提供了T集群账户。关于T集群的使用，请大家参考 [Run jobs in a cluster \[Advanced\]](/teaching/getting-started/cluster.md)
* 完成大作业需要大家自己编写一些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)
```

![](/files/1bAFGIITdOoTUKlGZBKw)

* 有了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表达量的例子。

![](/files/vY0gW5b0Kjc9MXvl146m)

### 对特征进行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 and samtools](/teaching/part-iii.-ngs-data-analyses/1.mapping/1.2-bedtools-samtools.md)


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://book.ncrnalab.org/teaching/part-v.-assignments/1.quiz_exrna/help.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.