栏目分类
热点资讯
你的位置:开发一个教育物联网软件需要多少钱 > 物联网软件开发公司 > 开发一个教育物联网软件需要多少钱 使用Seurat的v5来读取多个不是10x模范文献的单细胞形势

物联网软件开发公司

开发一个教育物联网软件需要多少钱 使用Seurat的v5来读取多个不是10x模范文献的单细胞形势

发布日期:2024-11-06 05:16    点击次数:110

前边咱们在 初试Seurat的V5版块 的推文内部演示了10x单细胞样品的模范3文献的读取开发一个教育物联网软件需要多少钱,并且在使用Seurat的v5来读取多个10x的单细胞转录组矩阵 的推文内部演示了多个10x单细胞样品的模范3文献的读取。

然而留住来了一个悬念, 等于若是咱们的单细胞转录组并不是10x的模范3文献,而是tsv或者csv或者txt等文本文献抒发量矩阵信息,就有点羁系了。接下来咱们以2020的著作:《Single-Cell Transcriptome Analysis Reveals Dynamic Cell Populations and Differential Gene Expression Patterns in Control and Aneurysmal Human Aortic Tissue》例如诠释,它的数据集是 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE155468

不错看到,作家这个技能给出来的是:

GSM4704931_Con4.txt.gz 9.2 MbGSM4704932_Con6.txt.gz 3.0 MbGSM4704933_Con9.txt.gz 10.0 MbGSM4704934_TAA1.txt.gz 7.7 MbGSM4704935_TAA2.txt.gz 5.8 MbGSM4704936_TAA3.txt.gz 7.2 MbGSM4704937_TAA4.txt.gz 12.5 MbGSM4704938_TAA5.txt.gz 11.7 MbGSM4704939_TAA6.txt.gz 8.1 MbGSM4704940_TAA7.txt.gz 18.7 MbGSM4704941_TAA8.txt.gz 6.4 Mb

是11个单细胞转录组样品,8 patients with ATAA (4 women and 4 men) and 3 controls (2 women and 1 man). 每个样品都是一个沉寂的txt文本文献贮蓄着其抒发量矩阵信息。

值得醒想法是这个2020的数据集还被2023的著作援用了,感兴致的不错去望望,Genome-wide association study of thoracic aortic aneurysm and dissection in the Million Veteran Program. Nat Genet 2023 Jul;55(7):1106-1115. PMID: 37308786

app

前边提到了,若是是莫得样品的txt沉寂读取后,再merge的技能成为的Seurat对象内部的各个样品的抒发量矩阵的分开的,就会导致所有这个词的后头的身手都失败。而它每个样品并不是10x单细胞样品的模范3文献,是以没宗旨使用前边的战略。

第一种步调是把每个样品的矩阵对都

每个样品的txt仍然是沉寂的读取,代码如下所示:

dir='GSE155468_RAW/' samples=list.files( dir ,pattern = 'gz')samples library(data.table)ctList = lapply(samples,function(pro){   # pro=samples[1]   print(pro)  ct=fread(file.path( dir ,pro),data.table = F)  ct[1:4,1:4]  rownames(ct)=ct[,1]  colnames(ct) = paste(gsub('.txt.gz','',pro),                       colnames(ct) ,sep = '_')  ct=ct[,-1]   return(ct)})

上头的代码复返了 ctList 这个list,它内部有每个单细胞样品的抒发量矩阵,然而每个样品的基因数目和细胞数目都是不一样的哦。然后提前把矩阵吞并之前需要最初把基因数目对都,吞并后才构建对象:

lapply(ctList, dim)tmp =table(unlist(lapply(ctList, rownames)))cg = names(tmp)[tmp==length(samples)]bigct = do.call(cbind,                lapply(ctList,function(ct){                   ct = ct[cg,]                   return(ct)                }))sce.all=CreateSeuratObject(counts =  bigct,                        min.cells = 5,                       min.features = 300)sce.allas.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])head(sce.all@meta.data, 10)table(sce.all@meta.data$orig.ident) 

不错看到,我这个技能作念了一个处分,等于每个样品的基因数目都对都了,若是某个基因在某个样品内部独有其实它不会被计议,因为我计议的是绝大部分基因。

因为多个样品吞并成为了一个超等大的抒发量矩阵,等于 bigct 这个变量,是以后头径直针对它来使用CreateSeuratObject函数去构建Seurat对象,等于无缺的下流分析的输入数据啦。

第二种步调是把矩阵归附成为10x的3文献

前边咱们指出来了,它每个样品并不是10x单细胞样品的模范3文献,每个样品都是一个沉寂的txt文本文献贮蓄着其抒发量矩阵信息,是以没宗旨使用前边的战略。然而,物联网app开发咱们其实不错说明这个txt文献去把它归附成为10x的3文献,早在2020-03-16其实我就写个一个简便的札记:抒发矩阵逆转为10X的模范输出3个文献,然而阿谁技能的代码稍许有点羁系,咱们其实不错把它写成一个函数,接下来让咱们演示一下吧。

每个样品的txt仍然是沉寂的读取,代码如下所示:

dir='GSE155468_RAW/' samples=list.files( dir ,pattern = 'gz')samples library(data.table)ctList = lapply(samples,function(pro){   # pro=samples[1]   print(pro)  ct=fread(file.path( dir ,pro),data.table = F)  ct[1:4,1:4]  rownames(ct)=ct[,1]  colnames(ct) = paste(gsub('.txt.gz','',pro),                       colnames(ct) ,sep = '_')  ct=ct[,-1]   return(ct)})

上头的代码复返了 ctList 这个list,它内部有每个单细胞样品的抒发量矩阵,然而每个样品的基因数目和细胞数目都是不一样的哦。接下来咱们构造一个自界说函数,把抒发量矩阵转为10x的3个文献,如下所示:

to10x <- function(ct)  {  write.table(data.frame(rownames(ct),rownames(ct)),file = 'features.tsv',              quote = F,sep = '\t',              col.names = F,row.names = F)  write.table(colnames(ct),file = 'barcodes.tsv',quote = F,              col.names = F,row.names = F)  file="matrix.mtx"  sink(file)  cat("%%MatrixMarket matrix coordinate integer general\n")  cat("%\n")  cat(paste(nrow(ct),ncol(ct),sum(ct>0),"\n"))   sink()  tmp=ct[1:5,1:4]  tmp  tmp=do.call(rbind,lapply(1:ncol(ct),function(i){    return(data.frame(row=1:nrow(ct),                      col=i,                      exp=ct[,i]))  }) )  tmp=tmp[tmp$exp>0,]  head(tmp)  write.table(tmp,file = 'matrix.mtx',quote = F,              col.names = F,row.names = F,append = T )}

比拟简便,接下来就针对前边的抒发量列表去轮回使用这个函数即可,如下所示:

 lapply(samples,function(pro){   # pro=samples[1]   pro=gsub('.txt.gz','',pro)  print(pro)  ct = ctList[[1]]  dir.create(pro)  setwd(pro)  to10x(ct)  setwd('../')  })

说真话,函数初始恶果如实有点低,不外不勤勉,反恰是进修的代码,咱们粗豪是仍是会遴选前边的矩阵吞并的模式,并不需要把抒发量矩阵转为10x的3个文献。成效后不错看到如下所示的文献夹结构:

│   ├── [ 160]  GSM4704935_TAA2│   │   ├── [115K]  barcodes.tsv│   │   ├── [291K]  features.tsv│   │   └── [ 95M]  matrix.mtx│   ├── [ 160]  GSM4704936_TAA3│   │   ├── [115K]  barcodes.tsv│   │   ├── [291K]  features.tsv│   │   └── [ 95M]  matrix.mtx

值得醒想法是每个样品这个技能内部的3文献其实是并莫得压缩,是以很滥用空间哦。并且因为这个技能我给出来的名字是features.tsv是以若是念念使用Seurat的Read10X读取,就需要把每个样品文献夹内部的3文献gz压缩一下哦!然后把每个样品的文献夹归纳整理到 outputs 文献夹内部,就不错使用如下所示的代码啦。

龙头分析:上期龙头开出1区号码05,与前期比对点位上升,本期关注龙头点位下降,推荐03。

library(Seurat)tmp = list.dirs('outputs')[-1]tmpct = Read10X(tmp) sce.all=CreateSeuratObject(counts = ct  ,                            min.cells = 5,                           min.features = 300,) 

如下所示的文献夹架构哦:

图片

文献夹架构

雷同的,只需有了sce.all对象开发一个教育物联网软件需要多少钱,后头的降维聚类分群等于咱们之前的代码即可啦。未来中午直播一下这个一说念的历程哈!

本站仅提供存储处事,所有这个词实际均由用户发布,如发现存害或侵权实际,请点击举报。