网站如何申请微信支付,官方网站welcome怎么注册,北京seo如何排名,wordpress改字号写在前边
虽然现在是高通量测序的时代#xff0c;但是GEO、ArrayExpress等数据库储存并公开大量的基因表达芯片数据#xff0c;还是会有大量的需求去处理芯片数据#xff0c;并且建模或验证自己所研究基因的表达情况#xff0c;芯片数据的处理也可能是大部分刚学生信的道友…写在前边
虽然现在是高通量测序的时代但是GEO、ArrayExpress等数据库储存并公开大量的基因表达芯片数据还是会有大量的需求去处理芯片数据并且建模或验证自己所研究基因的表达情况芯片数据的处理也可能是大部分刚学生信的道友入门R语言数据处理的第一次实战因此准备更新100个基因表达芯片或转录组高通量数据的处理。
数据信息检索
可以看到GSE126848是转录组高通量测序数据因此可以使用GEOquery包下载数据临床信息并且手动下载表达矩阵并整理 使用GEOquery包下载数据
using(tidyverse, GEOquery, magrittr, data.table, AnnoProbe, clusterProfiler, org.Hs.eg.db, org.Mm.eg.db)
注using是我写的函数作用是一次性加载多个R包不用写双引号并且不在屏幕上打印包的加载信息可以参考之前的推文using的定义函数名字using是在模仿Julia语言中的包加载函数
geo_accession - GSE126848
gset - GEOquery::getGEO(geo_accession, destdir ./, AnnotGPL F, getGPL F)
eSet - gset[[1]]
gpl - eSetannotation处理表型数据
这部分是很关键的可以筛选一下分组表型信息只保留自己需要的样本在这里只保留disease:ch1中healthy和NASH的样本作为后续分析的样本根据自己的研究目的筛选符合要求的样本
pdata - pData(eSet)
geo_accessiondescriptiondisease:ch1gender:ch1tissue:ch1GSM36152932683NAFLDMaleLiverGSM36152942685NAFLDMaleLiverGSM36152952687NAFLDMaleLiverGSM36152962689NAFLDFemaleLiverGSM36152972691NAFLDFemaleLiverGSM36152982693NAFLDMaleLiver
pdata %%dplyr::mutate(Sample geo_accession,Group case_when(diagnosis:ch1 HC ~ Control, diagnosis:ch1 NASH ~ Case, TRUE ~ NA),Age age (y):ch1,Sex str_to_title(gender:ch1),Stage fibrosis (stage):ch1) %%dplyr::filter(!is.na(Group)) %%dplyr::select(Sample, Group, Age, Sex)
fwrite(pdata, file str_glue({geo_accession}_pdata.csv))
处理表达谱数据
原始数据为Count值需要标准化为TPM并且基因名是Ensembl ID转换为Symbol基因名可以使用到我自己写的几个函数genekit、bioquest有需要可以联系我的公众号恩喜玛生物加入交流群
import pandas as pd
import genekit as gk
import bioquest as bqfdata pd.read_csv(GSE126848_Gene_counts_raw.txt.gz,sep\t,index_col0)
pdata pd.read_csv(GSE126848_pdata.csv,index_col0)
pdata.drop(columns[Sample2]).to_csv(GSE126848_pdata.csv)
fdata与pdata样本名统一这里使用了Python的字符串格式化方法
fdata fdata.loc[:,[{0:04}.format(x) for x in pdata.Sample2]]
fdata.columns pdata.index.to_list()
保存一份原始Count数据信息
fdata.to_csv(GSE126848_count.csv.gz)Count 转 TPM
fdata gk.countto(fdata, towhattpm, geneidEnsembl, speciesHuman)Ensembl ID转换为Symbol基因名
fdatagk.geneIDconverter(framefdata,from_idEnsembl,to_idSymbol,keep_fromFalse,gene_typeFalse,)去重复
根据每个基因表达量的中位数去除重复的基因
fdatabq.tl.unique_exprs(fdata)保存TPM基因表达量数据
fdata.to_csv(GSE126848_tpm.csv.gz)
