生信数据分析|scRNA-seq数据获取与转化
生信数据下载实在是太麻烦了!!
由于生信数据下载实在是太麻烦了!!网上的方法太多太冗杂,特写篇博客记录一下最近半年用到的方法。值得注意的是,以下部分代码不是原创(不过整理是原创的!),但是出处已经找不到了。此篇博客仅供学习交流作用。
1.单细胞数据获取
1.1借助包来下载GEO官网的数据
1.1.1方法一:利用GEOquery
以GSE213252为例,可以用GEOquery来做到。
if(!requireNamespace('BiocManager',quietly = TRUE))
+ install.packages('BiocManager')
BiocManager::install('GEOquery')
library(GEOquery)
cancer <- getGEO('GSE213252')
cancer1 <- cancer[[1]]
1.1.2方法二:利用GEOquery
以GSE148673为例,如果方法一行不通,可以试试运行下方代码。
suppressMessages(library(GEOquery))
library(Biobase)
options( 'download.file.method.GEOquery' = 'libcurl' )
gds <- getGEO("GSE148673" , destdir = "./",AnnotGPL = FALSE,getGPL = FALSE)
gds1 <- gds[[1]]
1.1.3方法三:更改国内镜像
#the function used in this method require "a file"
#"Your GSE may not be expression by array, or even not a GSE"
install.packages("devtools")
devtools::install_git("https://gitee.com/jmzeng/GEOmirror")
remotes::install_github("jmzeng1314/AnnoProbe")
library(GEOquery)
library(GEOmirror)
library(AnnoProbe)
eSet=AnnoProbe::geoChina('GSE131928')
eSet
eSet=eSet[[1]]
phenoDat <- pData(eSet)
1.2手动下载并导入数据
顾名思义,先去GEO官网下载对应数据,然后再导入r或python里。以下根据不同的数据格式提供多种方法参考。
1.2.1方法一:features\barcodes\matrix
无论是r还是python都要记得把这三个文件都改名,并放在同一个文件夹下。如下图所示。
【r版导入成Seurat对象】
library(Seurat)
setwd("/Users/mhuang/Downloads")
dir <- "./data/"
list.files(dir)
counts <- Read10X(data.dir = dir)
class(counts)
scRNA <- CreateSeuratObject(counts = counts)
scRNA
colnames(scRNA@meta.data$cell_types) #按理说是放在这里的
【python版导入成adata】
python直接用scanpy包的read_10x_mtx函数导入即可。以我上方的文件夹截图为例,由于我是control_rep1文件夹里存放着这三个文件,因此路径只要写到control_rep1就好。
import scanty as sc
adata = sc.read_10x_mtx("/Users/mhuang/Downloads/control_rep1")
如果python版导入的时候出现问题了,那可能是文件本身不是标准格式,是需要更改的。以我上方的文件夹截图为例,由于文件太大, 因此都是以gz结尾的压缩包。如果我导入失败,那么我首先应该在终端上解压文件(以matrix文件为例,解压方式如下图所示),然后检查它是否符合标准格式,若不符合,就更改。
gunzip -c matrix.mtx.gz | head
标准格式可以查看下方的参考文献:
seurat读取文件的格式 10x文件内容 mtx格式scanpy.read_10x_mtx scanpy读取10x格式_10x单细胞mtx文件中的基因列和细胞列互换了,如何换回来-CSDN博客
1.2.2方法二:csv
方法二和方法三单纯是我记不住函数,所以全部都写出来一下!
xb=read.csv("/Users/mhuang/Downloads/counts.table.csv")
#两个参数的设置和read.table函数一样
1.2.3方法三:txt\tsv
file=read.table("/Users/mhuang/Downloads/counts.table.tsv",fill = TRUE,sep="\t")
#如果有缺失值就可以加上这个指标
colnames(file)=file[1,] #将第一行作为列名(或者read.table加上参数header=TRUE)
rownames(file)=file[,1] #将第一列作为行名(或者read.csv加上参数row.names=1)
file <- file[-1,] #删除第一行
1.2.4方法四:h5\hdf5
这边仍然提供两种方法,如下所示。
#1.
BiocManager::install("rhdf5")
library(rhdf5)
h5_file= H5Fopen("new.h5")
####new.h5文件内创建了一个组(group1_mat),组内又创建了df和matrix两个层级用以保存矩阵和数据框
h5dump(h5_file,load=FALSE)
#2.
BiocManager::install("rhdf5")
library(rhdf5)
library(Seurat)
data_sample <- Read10X_h5("xxx.h5")
data_seurat <- CreateSeuratObject(data_sample,project = "data_sample")
2.单细胞数据转化
2.1h5ad转为Seurat
2.1.1方法一:利用seurat-disk包来转化
这个方法是最常见的教程,但是我基本上没运行成功过。。
remotes::install_github("mojaveazure/seurat-disk")
library(SeuratDisk)
Convert("xxxxx.h5ad", dest="h5seurat",
assay = "RNA",
overwrite=F)
seurat_object <- LoadH5Seurat("xxxxx.h5seurat")
2.1.2方法二:先转为h5和tsv文件,再导入r中
这个方法百试百灵,应该还能运用到st数据(空转数据)中,不过我没试过,有需要可以试试。主要是分为两步:
第一步,把注释信息和表达矩阵在python中分离出来。
import h5py
import pandas as pd
import scanpy as sc
adata = sc.read("/file_root/PBMC_Batch_test.h5ad")
data_R=pd.DataFrame(data=adata.X, #本来要加个todense()的,但是这边不知道为啥已经是密集矩阵了
index=adata.obs_names,
columns=adata.var_names)
data_R.to_hdf("/file_root/process_test_res.h5","data_R")
meta_data=pd.DataFrame(data=adata.obs)
meta_data.to_csv('/file_root/process_test_meta.tsv',sep="\t")
第二步,导入进r并创造seurat对象。
data_R <- h5read("/file_root/process_test_res.h5","data_R")
mat <- data_R$block0_values
rownames(mat) <- data_R$axis0
colnames(mat) <- data_R$axis1
mat <- Matrix(mat, sparse = TRUE) #转为稀疏矩阵
meta_data <- read.table('/file_root/process_test_meta.tsv',sep="\t",header=T,row.names=1)
facs_test <- CreateSeuratObject(mat,project='Spatial',meta.data=meta_data)
2.2SingleCellExperiment对象转为h5ad
主要是用python实现的转化,其中利用rpy2包来运行r代码。
import numpy as np
import pandas as pd
import scanpy as sc
import anndata as ad
from rpy2.robjects.packages import importr
importr("Seurat")
r('''
sce = readRDS('dataset_C(1).Rds')
''')
obs = pd.DataFrame()
celltype = r('sce@colData@listData[["X"]]')
obs['cell.type']= pd.DataFrame(celltype)
X = r('sce@assays@data@listData[["counts"]]')
X = np.array(X).T
index = r('sce@colData@rownames')
adata = ad.AnnData(X, obs=obs, dtype='int32')
adata.obs.index = index
adata.write_h5ad('adata4.h5ad')
2.3 SingleCellExperiment对象转为Seurat
library(Seurat)
colnames(pbmc.sce) <- make.unique(colnames(pbmc.sce))#要确保列名唯一
pbmc.seurat <- as.Seurat(pbmc.sce)
2.4 Seurat转为h5ad
其实是2.1.1节的函数反过来运行即可。
library(SeuratDisk)
library(Seurat)
SaveH5Seurat(seurat,filename = "dlpfc.h5Seurat")
Convert("dlpfc.h5Seurat",dest="h5ad")
更多推荐
所有评论(0)