python 改行学it 生信数据分析｜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')

相关文章