recount3超大规模可用转录组数据集

2023-12-13

随着测序数据的积累，如何复用这些数据是一个挑战。recount项目，目前是recount3，共收集了8,679人和10,088小鼠的数据集，超过750,000个样本，经过统一处理（uniformly processed），得到gene或exon的表达以及exon-exon junction的数据。好多年前，我就了解过recount项目，很奇怪很少有介绍的。

一，recount对所有属于进行了uniformly processed，避免了分析流程的bias；

二，recount提供了超大规模的预处理之后的数据，直接拿来用，避免研究人员从原始数据分析；

三，recount提供了简单易用的工具，方便研究人员下载和处理数据。

方法1：下载TCGA-OV为例，检索过滤然后下载

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26


library(recount3)

# 同 https://jhubiostatistics.shinyapps.io/recount3-study-explorer/
# 可以看到project_home和project，包括TCGA,GTEX和SRA
human_projects <- available_projects()

proj_info <- subset(
  human_projects,
  project == "OV" & project_type == "data_sources"
)

rse_gene_tcga_ov <- create_rse(proj_info)

#counts data
assays(rse_gene_tcga_ov)$counts <- transform_counts(rse_gene_tcga_ov)
# ## Compute TPMs
assays(rse_gene_tcga_ov)$TPM <- recount::getTPM(rse_gene_tcga_ov, length_var = "score")
# ## Check TPM. Should all be equal to 1
colSums(assay(rse_gene_tcga_ov, "TPM")) / 1e6 


# View sample annotation
View(data.frame(colData(rse_gene_tcga_ov)))

# View gene annotation
View(data.frame(exp$tcga.ov.expr@rowRanges))

方法2：直接选中数据集，生成下载code

在这个网站选中想要下载的数据集，https://jhubiostatistics.shinyapps.io/recount3-study-explorer/，网站下方会显示下载的code。注释不一定是26，还可以是RefSeq v109，Gencode v29等。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14


rse_gene_tcga_ov = recount3::create_rse_manual(
    project = "OV",
    project_home = "data_sources/tcga",
    organism = "human",
    annotation = "gencode_v26",
    type = "gene"
)

#counts data
assays(rse_gene_tcga_ov)$counts <- transform_counts(rse_gene_tcga_ov)
# ## Compute TPMs
assays(rse_gene_tcga_ov)$TPM <- recount::getTPM(rse_gene_tcga_ov, length_var = "score")
# ## Check TPM. Should all be equal to 1
colSums(assay(rse_gene_tcga_ov, "TPM")) / 1e6

Kingfisher下载SRA数据

2023-11-23

Default Category

我知道一个SRP的编号，里面有我想要下载的数据，我想根据SRP编号快速下载数据，查到了Kingfisher这个工具。

https://github.com/wwood/kingfisher-download

文档：https://wwood.github.io/kingfisher-download/

安装：pip install kingfisher

主要有三个模块，get、extract、annotate

get

1
2
3
4
5
6


kingfisher get -r ERR1739691 -m ena-ascp
# 可以指定列表--run-identifiers-list
# 指定project编号SRP，--bioprojects
# 指定下载方法，--download-methods，可以指定多种，程序会一个方法一个方法的试，包括ena-ascp，ena-ftp，prefetch，aws-http，aws-cp，gcp-cp
# 指定线程数目，--download-threads
# 指定ascp需要的key路径，--ascp-ssh-key

比对、定量多种类型小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/"))