生信小白开始接活——文本处理（二）

上篇处理了最简单的单细胞公共数据，即1个样本对应1个矩阵文件（tsv），读取统计就行了，但是GEO公共数据库中还存在如下的情况：例如GSE178318，15个样本的信息全部混在一块了，针对这种情况，我们怎么统计指定样本的线粒体基因分布，中值基因，中值UMI这些指标，逻辑就是拆分样本，如何拆分，下面会娓娓道来~~~

1.png

打开GEO界面，红框处6个样本的数据需要统计：

2.png

15个样本的数据见下矩阵：

3.png

拆分样本

1、下载3个矩阵文件

wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE178nnn/GSE178318/suppl/GSE178318_barcodes.tsv.gz
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE178nnn/GSE178318/suppl/GSE178318_genes.tsv.gz
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE178nnn/GSE178318/suppl/GSE178318_matrix.mtx.gz
barcodes.tsv.gz   features.tsv.gz  matrix.mtx.gz #将下载下来的3个文件重新命名，命令为mv 原名 新名，方便直接用Read10X函数读取

2、在 Seurat 中，对于 Cellranger 数据的导入，它提供了两个函数 Read10X 和 Read10X_h5。前者用于读取 Cellranger 生成的 3 个文件：barcodes.tsv，features.tsv，matrix.mtx，后者则用于读取生成的 xxxxxx.h5文件，这里我们应用前面一个函数就行

R #linux系统调用R
library(Seurat) #加载seurat软件
a <- Read10X("~/SingleronTest/wuxuan/leo") #读取3个文件
a <- CreateSeuratObject(a) #创建Seurat对象
head(a@meta.data) #查看读取的数据前6行，结果如下展示
 orig.ident nCount_RNA nFeature_RNA
AAACCTGAGAAACCTA_COL07_CRC SeuratProject       5565         1241
AAACCTGAGACAATAC_COL07_CRC SeuratProject       5892         1845
AAACCTGAGACGCAAC_COL07_CRC SeuratProject       3372          933
AAACCTGAGCAGATCG_COL07_CRC SeuratProject        945          548
AAACCTGAGCTATGCT_COL07_CRC SeuratProject       1565          571
AAACCTGAGGACATTA_COL07_CRC SeuratProject       1865          816
head(a@meta.data, 8) #查看前8行的命令

3、正式拆分，使用tidyverse函数

library(tidyverse) #加载tidyverse

4、拆分脚本

rownames(a@meta.data) -> a@meta.data$rowLeo #反向赋值，增加了1列
head(a@meta.data) #展示如下
                              orig.ident nCount_RNA nFeature_RNA
AAACCTGAGAAACCTA_COL07_CRC SeuratProject       5565         1241
AAACCTGAGACAATAC_COL07_CRC SeuratProject       5892         1845
AAACCTGAGACGCAAC_COL07_CRC SeuratProject       3372          933
AAACCTGAGCAGATCG_COL07_CRC SeuratProject        945          548
AAACCTGAGCTATGCT_COL07_CRC SeuratProject       1565          571
AAACCTGAGGACATTA_COL07_CRC SeuratProject       1865          816
                                               rowLeo
AAACCTGAGAAACCTA_COL07_CRC AAACCTGAGAAACCTA_COL07_CRC
AAACCTGAGACAATAC_COL07_CRC AAACCTGAGACAATAC_COL07_CRC
AAACCTGAGACGCAAC_COL07_CRC AAACCTGAGACGCAAC_COL07_CRC
AAACCTGAGCAGATCG_COL07_CRC AAACCTGAGCAGATCG_COL07_CRC
AAACCTGAGCTATGCT_COL07_CRC AAACCTGAGCTATGCT_COL07_CRC
AAACCTGAGGACATTA_COL07_CRC AAACCTGAGGACATTA_COL07_CRC
str_split(a$rowLeo,"_")  #拆的结果如下展示
[[99997]]
[1] "CCGTACTAGGGTATCG" "COL17"            "LM"              

[[99998]]
[1] "CCGTACTCAAGCCGCT" "COL17"            "LM"              

[[99999]]
[1] "CCGTACTCAAGGGTCA" "COL17"            "LM"   

unlist(str_split(a$rowLeo,"_")) #拆的结果如下展示 
[99988] "AACTCAGCAAGGTTTC" "COL12"            "CRC"             
[99991] "AACTCAGCACTCAGGC" "COL12"            "CRC"             
[99994] "AACTCAGGTATAGTAG" "COL12"            "CRC"             
[99997] "AACTCAGGTCATATCG" "COL12"            "CRC"             
 [ reached getOption("max.print") -- omitted 320844 entries ]
head(str_split(a$rowLeo,"_",simplify=T)) #展示前6行
     [,1]               [,2]    [,3] 
[1,] "AAACCTGAGAAACCTA" "COL07" "CRC"
[2,] "AAACCTGAGACAATAC" "COL07" "CRC"
[3,] "AAACCTGAGACGCAAC" "COL07" "CRC"
[4,] "AAACCTGAGCAGATCG" "COL07" "CRC"
[5,] "AAACCTGAGCTATGCT" "COL07" "CRC"
[6,] "AAACCTGAGGACATTA" "COL07" "CRC"
head(str_split(a$rowLeo,"_",simplify=T) [,2]) #展示[,2]的前6行
[1] "COL07" "COL07" "COL07" "COL07" "COL07" "COL07"
str_split(a$rowLeo,"_",simplify=T) [,2] -> a@meta.data$Sample #反向赋值
str_split(a$rowLeo,"_",simplify=T) [,3] -> a@meta.data$Tissue #反向赋值
head(a@meta.data)
                              orig.ident nCount_RNA nFeature_RNA
AAACCTGAGAAACCTA_COL07_CRC SeuratProject       5565         1241
AAACCTGAGACAATAC_COL07_CRC SeuratProject       5892         1845
AAACCTGAGACGCAAC_COL07_CRC SeuratProject       3372          933
AAACCTGAGCAGATCG_COL07_CRC SeuratProject        945          548
AAACCTGAGCTATGCT_COL07_CRC SeuratProject       1565          571
AAACCTGAGGACATTA_COL07_CRC SeuratProject       1865          816
                                               rowLeo Sample Tissue
AAACCTGAGAAACCTA_COL07_CRC AAACCTGAGAAACCTA_COL07_CRC  COL07    CRC
AAACCTGAGACAATAC_COL07_CRC AAACCTGAGACAATAC_COL07_CRC  COL07    CRC
AAACCTGAGACGCAAC_COL07_CRC AAACCTGAGACGCAAC_COL07_CRC  COL07    CRC
AAACCTGAGCAGATCG_COL07_CRC AAACCTGAGCAGATCG_COL07_CRC  COL07    CRC
AAACCTGAGCTATGCT_COL07_CRC AAACCTGAGCTATGCT_COL07_CRC  COL07    CRC
AAACCTGAGGACATTA_COL07_CRC AAACCTGAGGACATTA_COL07_CRC  COL07    CRC
tail(a@meta.data) 
                               orig.ident nCount_RNA nFeature_RNA
TTTGTCAGTTCCACGG_COL18_PBMC SeuratProject       5578         1346
TTTGTCAGTTGGGACA_COL18_PBMC SeuratProject       3485         1192
TTTGTCATCATGTGGT_COL18_PBMC SeuratProject       4894         1445
TTTGTCATCCGCATCT_COL18_PBMC SeuratProject       4051         1111
TTTGTCATCGGAATCT_COL18_PBMC SeuratProject       2617         1018
TTTGTCATCGGTTAAC_COL18_PBMC SeuratProject       2976         1235
                                                 rowLeo Sample Tissue
TTTGTCAGTTCCACGG_COL18_PBMC TTTGTCAGTTCCACGG_COL18_PBMC  COL18   PBMC
TTTGTCAGTTGGGACA_COL18_PBMC TTTGTCAGTTGGGACA_COL18_PBMC  COL18   PBMC
TTTGTCATCATGTGGT_COL18_PBMC TTTGTCATCATGTGGT_COL18_PBMC  COL18   PBMC
TTTGTCATCCGCATCT_COL18_PBMC TTTGTCATCCGCATCT_COL18_PBMC  COL18   PBMC
TTTGTCATCGGAATCT_COL18_PBMC TTTGTCATCGGAATCT_COL18_PBMC  COL18   PBMC
TTTGTCATCGGTTAAC_COL18_PBMC TTTGTCATCGGTTAAC_COL18_PBMC  COL18   PBMC

5、按样本（Sample）提取

Leo <- SplitObject(a,split.by='Sample') 
Leo
$COL07
An object of class Seurat 
33694 features across 33100 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

$COL12
An object of class Seurat 
33694 features across 26310 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

$COL15
An object of class Seurat 
33694 features across 13950 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

$COL16
An object of class Seurat 
33694 features across 14200 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

$COL17
An object of class Seurat 
33694 features across 25735 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

$COL18
An object of class Seurat 
33694 features across 26986 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

6、按组织（Tissue）提取

Xuan <- SplitObject(a,split.by='Tissue') 
Xuan
$CRC
An object of class Seurat 
33694 features across 55735 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

$LM
An object of class Seurat 
33694 features across 57596 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

$PBMC
An object of class Seurat 
33694 features across 26950 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

6、因为只要属于CRC的6个样本，展示下文本结果便于理解

head(Xuan$CRC@meta.data)
                              orig.ident nCount_RNA nFeature_RNA
AAACCTGAGAAACCTA_COL07_CRC SeuratProject       5565         1241
AAACCTGAGACAATAC_COL07_CRC SeuratProject       5892         1845
AAACCTGAGACGCAAC_COL07_CRC SeuratProject       3372          933
AAACCTGAGCAGATCG_COL07_CRC SeuratProject        945          548
AAACCTGAGCTATGCT_COL07_CRC SeuratProject       1565          571
AAACCTGAGGACATTA_COL07_CRC SeuratProject       1865          816
                                               rowLeo Sample Tissue
AAACCTGAGAAACCTA_COL07_CRC AAACCTGAGAAACCTA_COL07_CRC  COL07    CRC
AAACCTGAGACAATAC_COL07_CRC AAACCTGAGACAATAC_COL07_CRC  COL07    CRC
AAACCTGAGACGCAAC_COL07_CRC AAACCTGAGACGCAAC_COL07_CRC  COL07    CRC
AAACCTGAGCAGATCG_COL07_CRC AAACCTGAGCAGATCG_COL07_CRC  COL07    CRC
AAACCTGAGCTATGCT_COL07_CRC AAACCTGAGCTATGCT_COL07_CRC  COL07    CRC
AAACCTGAGGACATTA_COL07_CRC AAACCTGAGGACATTA_COL07_CRC  COL07    CRC
head(Xuan$CRC@meta.data$Sample,100000)  #查看CRC里面的Sample
[55693] "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18"
[55702] "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18"
[55711] "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18"
[55720] "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18"
[55729] "COL18" "COL18" "COL18" "COL18" "COL18" "COL18" "COL18"
table(Xuan$CRC@meta.data$Sample)
COL07 COL12 COL15 COL16 COL17 COL18 
17000 10750  4250  6000  9035  8700 
head(Xuan$CRC@meta.data$Sample,17001)
[16966] "COL07" "COL07" "COL07" "COL07" "COL07" "COL07" "COL07" "COL07" "COL07"
[16975] "COL07" "COL07" "COL07" "COL07" "COL07" "COL07" "COL07" "COL07" "COL07"
[16984] "COL07" "COL07" "COL07" "COL07" "COL07" "COL07" "COL07" "COL07" "COL07"
[16993] "COL07" "COL07" "COL07" "COL07" "COL07" "COL07" "COL07" "COL07" "COL12"

7、第二次拆分

Xuan2 <- SplitObject(Xuan$CRC,split.by='Sample') 
Xuan2
$COL07
An object of class Seurat 
33694 features across 17000 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

$COL12
An object of class Seurat 
33694 features across 10750 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

$COL15
An object of class Seurat 
33694 features across 4250 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

$COL16
An object of class Seurat 
33694 features across 6000 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

$COL17
An object of class Seurat 
33694 features across 9035 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

$COL18
An object of class Seurat 
33694 features across 8700 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

统计数据

至此拆分结束，接下来就到了统计线粒体基因分布，中值基因，中值UMI这些指标的时候了

Xuan2$COL12 -> COL12 #反向赋值
median(COL12 @meta.data$nFeature_RNA)
[1] 804
median(COL12@meta.data$nCount_RNA)
[1] 2180
COL12[["percent.mt"]] <- PercentageFeatureSet( COL12, pattern = "^MT-")
median(COL12@meta.data$percent.mt)
[1] 3.362562
MT5 <- subset(COL12,percent.mt < 5)
MT5
An object of class Seurat 
33694 features across 8608 samples within 1 assay 
Active assay: RNA (33694 features, 0 variable features)

OK，大工告成，你学废了吗？

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 212,185评论 6赞 493
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 90,445评论 3赞 385
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 157,684评论 0赞 348
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 56,564评论 1赞 284
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 65,681评论 6赞 386
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 49,874评论 1赞 290
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 39,025评论 3赞 408
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 37,761评论 0赞 268
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 44,217评论 1赞 303
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 36,545评论 2赞 327
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 38,694评论 1赞 341
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 34,351评论 4赞 332
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 39,988评论 3赞 315
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 30,778评论 0赞 21
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,007评论 1赞 266
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 46,427评论 2赞 360
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 43,580评论 2赞 349

生信小白开始接活——文本处理（二）

拆分样本

统计数据

推荐阅读更多精彩内容