本节中我们提供一些和大作业的编程实现相关的补充材料。
0) 编程环境
完成大作业需要大家自己编写一些python,R和shell脚本。我们推荐大家根据目的的不同选择合适的编程环境。
有些脚本需要一定的计算量,并且会反复使用。这类脚本一经确定,一般就不会反复改动,而是会通过更换命令行参数改变输入输出以及脚本的行为。在实际计算中,运行这类脚本就不太适合使用交互式的工具,建议大家在进行实际计算前做好测试,在计算时用好python和R的命令行参数。
画图本身一般对计算资源要求较低,并且很可能需要反复更换不同的展示方法,尝试不同的参数等,所以在本地用jupyter notebook或rstudio这样的交互式的IDE会比较方便。建议大家把在服务器上计算好的,可以用来画图的结果以文本文档的形式保存,从而将有一定计算量的脚本和画图的脚本独立开来。
有一些IDE,vscode和pycharm,提供了远程连接的功能。如果有同学习惯使用IDE,又想利用远程的计算资源,可以尝试这些工具。
在T集群上,我们已经装好了一些常用工具,位于/data/2022-bioinfo-shared/softwares
。大家可以将需要的工具添加到自己的环境变量中。
我们提供了预装了一些常用package的python和R的miniconda环境,可以通过如下方式进入:
Copy 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
两个子目录下。在使用中大家不需要 把这些数据都复制到自己的目录下,为了方便在自己的目录下访问,只需要建立一个符号链接即可:
Copy 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
Copy # 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
):
Copy 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中所占的比例。大家可以参考我们提供的脚本:
Copy # 按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:
Copy 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。大家也可以参考我们提供的脚本来进行数据整理:
Copy 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的例子,鼓励大家多做其他尝试:
Copy 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等。
Copy y <- edgeR::calcNormFactors(y,method="TMM")
cpm.matrix <- edgeR::cpm(y)
Copy library(EDASeq)
plotRLE(cpm.matrix)
去除批次效应
有很多R package都提供了去除批次效应的方法。以下提供了两个例子:
Copy 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)
Copy 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
4) 特征选择
请大家注意,特征选择原则上和模型拟合一样,不能用到测试集的label。
5) 模型评估和特征解释
6) peak calling
对于除了miRNA和piRNA之外的其他RNA类型,对于每个样本的每个bam文件,我们可以先用samtools去除duplication
Copy #!/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文件:
Copy #以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和它依赖的动态库加入环境变量:
Copy 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
Copy #以lncRNA为例
Piranha -b 20 -s lncRNA.tags.bed > lncRNA.peaks.bed
我们把bam文件和计算出的peak load到IGV上进行可视化,进而对peak的形状有一些直观的认识
根据peak calling的结果,我们再用去duplication之前的bam文件统计落在每一个peak内部的reads数量。以下是一个样本lncRNA的例子:
Copy 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的坐标上)或者自行定义一些区域
Last updated 8 months ago