比对、定量多种类型小RNA

2023-09-22
Default Category

前面讲了用工具定量小RNA，但要么软件安装困难，要么维护不好给种错误，所以自己也搭建（抄）了一套流程。主要思路是bowtie比对，HTSeq定量。其实也是利用了HTSeq定量需要gff文件的特点，比对到全基因组后只需要准备好gtf文件即可。

bowtie和bowtie的index文件我已经有了，是在miRDeep2流程中定量miRNA准备的，本文主要介绍定量其他种类的非编码小RNA。

比对

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

zcat input.fastq.gz | bowtie --seedlen 10 -p 48 -v2 -m20 --best -S --strata --chunkmbs 8000 $BOWTIE_IND - output.sam

-M 1           like -m, but reports 1 random hit (MAPQ=0); requires --best
- 为zcat的输入，bowtie不识别gz，所以用zcat管道流入bowtie
-p 48 线程
-v report end-to-end hits w/ <=v mismatches; ignore qualitie
-m20 丢弃超过20次multi-map的reads
--best --strata multi-map的reads只报告最好的比对
--chunkmbs 8000 best-first搜索的最大内存（M）
-S 输出SAM格式
BOWTIE_IND bowtie的索引

SAM文件处理

1

samtools sort --threads 10 out.sam -o out_sorted.bam --output-fmt BAM && samtools index out_sorted.bam

gtf文件准备

snRNA,snoRNA等（来自gencode的注释）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

axel -n 10 https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.annotation.gtf.gz

zcat gencode.v44.annotation.gtf.gz | grep "snoRNA" > snoRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "Mt_rRNA" > Mt_rRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "Mt_tRNA" > Mt_tRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "rRNA" > rRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "snRNA" > snRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "sRNA" > sRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "scaRNA" > scaRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "miRNA" > miRNA.gtf

tRNA

1
2
3

# 也来自gencode，tRNA genes predicted by ENSEMBL on the reference chromosomes using tRNAscan-SE
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.tRNAs.gtf.gz
zcat gencode.v44.tRNAs.gtf.gz > tRNA.gtf
阅读更多

COMPSRA定量多种类型小RNA

2023-09-21
Default Category

COMPSRA: a COMprehensive Platform for Small RNA-seq data AnalySis

软件和数据准备

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

# 下载，注意COMPSRA需要JRE
mkdir COMPSRA
cd COMPSRA
wget https://github.com/cougarlj/COMPSRA/raw/master/COMPSRA_V1.0.3.zip
unzip COMPSRA_V1.0.3.zip

# 手动安装star
cd bundle_v1/plug/star
wget -c wget -c https://github.com/alexdobin/STAR/archive/refs/tags/2.5.3a.zip
unzip 2.5.3a.zip


cd ../../..

# 下载注释数据，如果下载的文件不对，可以去https://github.com/cougarlj/COMPSRA上bundle1_V1/prebuilt_db中下载
java -jar COMPSRA.jar -tk -dr -ck miRNA_hg38,piRNA_hg38,tRNA_hg38,snoRNA_hg38,snRNA_hg38,circRNA_hg38

# 下载人基因组
java -jar COMPSRA.jar -tk -dr -ck star_hg38
阅读更多

安装R包sf、lwgeom等报错

2023-09-14
Default Category

安装R包sf、lwgeom的时候，报错

1
2

configure: error: proj_api.h not found in standard or given locations.
configure: error: libproj not found in standard or given locations.

要想办法让R在安装包的时候知道系统已经有了对应的library，我是先创建了一个sysR的conda环境，然后在虚拟环境下装对应的库，可能安多了-_-||

1
2
3
4
5
6
7

mamba install -c conda-forge gdal
mamba install -c conda-forge proj
mamba install -c conda-forge libgdal
mamba install -c r r-sf
mamba install -c conda-forge proj4
mamba install -c conda-forge r-lwgeom
mamba install -c conda-forge r-proj4

然后为了保险起见，我在bashrc中增加了LD的路径（我只是为了确保系统能找到library的路径，添件了感觉LD_LIBRARY_PATH，事后这个地方可以不用设置）。参见下面PS2，我git push遇到了问题，所以这里建议直接在Terminal里面export，不要写在bashrc中。

1
2

LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/master/zhu_zhong_xu/miniconda3/lib/:/master/zhu_zhong_xu/miniconda3/envs/sysR/lib
export LD_LIBRARY_PATH

在conda安装library的时候，会安装R base。但我用的是Rstudio server下的R，也就是系统级别的R，我为了在虚拟环境下调用系统的R，把conda/env/sysR/bin/R删掉了（切记），这样找不到环境下的R时，就会调用系统的R（当然可以用系统R的绝对路径）。然后在R里面安装对应的包，并指定库的文件夹。

1
2
3

install.packages("sf", configure.args=c("--with-proj-include=~/conda/envs/sysR/include, --with-proj-lib=~/conda/envs/sysR/lib/"))

install.packages("lwgeom", configure.args=c("--with-proj-include=~/conda/envs/sysR/include, --with-proj-lib=~/conda/envs/sysR/lib/"))
阅读更多

JC-ELMER一个基于甲基化和转录组进行调控分析的R包

2023-09-06
Default Category

**ELMER (Enhancer Linking by Methylation/Expression Relationships)**是基于甲基化数据和转录组数据，分析远端甲基化位点调控基因表达的R包，包里也整合了TCGA的数据，也可以自己构建MAE对象，快速进行分析，识别差异甲基化的位点，与基因表达显著相关的甲基化位点和相关区域的motif。

DNA 甲基化可用于识别肿瘤和其他原发性疾病组织中转录增强子和其他顺式调控模块（CRM，cis-regulatory modules）的功能变化。 R/Bioconductor 包 ELMER（通过甲基化/表达关系增强连接）提供了一种系统方法，通过结合来自同一组样本的甲基化和基因表达数据来重建基因调控网络 (GRN，gene regulatory networks)。 ELMER 使用 CRM 的甲基化变化作为网络调控的基础，利用相关性分析将它们与上游主调控 (MR，master regulator) 转录因子和下游目标基因相关联。

ELMER 分析有 5 个主要步骤：

  • 识别 HM450K 或 EPIC 阵列上的远端probe（集成了注释，很方便的拉取远端probe的信息）
  • 识别两组之间 DNA 甲基化水平显着不同的远端probe（常规的甲基化差异分析）
  • 识别差异甲基化远端探针的假定靶基因（通关相关性分析，关联probe和gene）
  • 识别远端探针的富集motif，这些基序具有显着差异甲基化并与假定的目标基因相关（识别motif）。
  • 识别其表达与富集motif处的 DNA 甲基化相关的调节性 TF转录因子

ELMER能做的分析，参考手册的plot页面，例如https://www.bioconductor.org/packages/release/bioc/vignettes/ELMER/inst/doc/plots_scatter.html

阅读更多

校正混杂因素

2023-08-17
Default Category

混杂因素亦称混杂因子或外来因素，指与研究因素和研究疾病均有关，若在比较的人群组中分布不均，可以歪曲（掩盖或夸大）因素与疾病之间真正联系的因素。

校正变量的方法很简单，只需要校正的变量和要分析的变量共同纳入方程即可，但是最好在纳入方程前对于自变量能有一个初筛即根据资料的特点和文献复习的情况，只纳入可能有关的,对于初筛p值特别大的最好不要纳入方程以免方程出现不稳定。

混杂因素的影响以及校正可以参考这个post：https://www.r-bloggers.com/2020/09/correcting-for-confounded-variables-with-glms/

阅读更多