Seurat로 cell type 결정하기
추천글 : 【생물정보학】 생물정보학 분석 목차
1. 정식 튜토리얼 [본문]
2. 커스토마이징 [본문]
3. 트러블슈팅 [본문]
a. Cell Type Classification Pipeline
※ 최근에 업그레이드 된 Seurat v5.0 이전의 문서입니다.
※ Seurat는 조르주 쇠라에서 따온 명칭입니다.
※ Seurat의 발음은 '스-라' (fr; 위) 혹은 '스-랫' (eng; 아래)에 가까움
1. 정식 튜토리얼 (R STUDIO) (ref1, ref2) [목차]
library(dplyr)
install.packages("Seurat")
library(Seurat)
### import input_data(cellranger count output과 동일) ###
pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")
### setup the seurat object ###
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
### Preprocessing(QC) ###
### Low quality cell, mitochondria genome percent check
pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
### Quality Control
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 15)
### normalizing the data ###
pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
### identification of highly variable features ###
pbmc <- FindVariableFeatures(object = pbmc, selection.method = "vst", nfeatures = 2000)
### scaling data ###
all.genes <- rownames(x = pbmc)
pbmc <- ScaleData(object = pbmc, features = all.genes)
### perform linear dimensional reduction ###
pbmc <- RunPCA(object = pbmc, features = VariableFeatures(object = pbmc))
### determine the "dimensionality" of the dataset
pbmc <- JackStraw(object = pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(object = pbmc, dims = 1:20)
### cluster the cells ###
pbmc <- FindNeighbors(object = pbmc, dims = 1:10)
pbmc <- FindClusters(object = pbmc, resolution = 0.5)
### Run non-linear dimensional reduction(tSNE) ###
pbmc <- RunUMAP(object = pbmc, dims = 1:10)
DimPlot(object = pbmc, reduction = "umap")
pbmc <- RunTSNE(object = pbmc, dims = 1:10)
DimPlot(object = pbmc, reduction = "tsne")
### Finding differentially expressed features(biomarkers) ###
markers <- FindAllMarkers(object = pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
### Maybe you need this kind of data ###
umap <- DimPlot(object = pbmc, reduction = "umap")
tsne <- DimPlot(object = pbmc, reduction = "tsne")
names(umap)
# [1] "data" "layers" "scales" "mapping" "theme"
# [6] "coordinates" "facet" "plot_env" "labels" "guides"
write.csv(markers, "markers.csv")
write.csv(umap$data, "umap_data.csv")
write.csv(tsne$data, "tsne_data.csv")
Figure. 1. DimPlot으로 나타낸 cell type에 따른 클러스터
⑴ 예제 파일
⑵ scRNA-seq data는 read count를 총 transcript로 나눈 뒤 10,000을 곱한 것을 나타냄
⑶ CreateSeuratObject
① min.features : 이 값이 작아야 누락되는 데이터가 적어짐
⑷ FindVariableFeatures 함수에서 nfeatures 파라미터를 통해 가장 발현율이 높은 2,000개의 유전자들이 선택됨
① nfeatures : 이 값이 커야 클러스터링이 잘 됨. 일반적으로 2,000 이상이면 클러스터링에 큰 영향을 주지 않음
⑸ NormalizeData
① scRNA-seq data는 log-normalization을 통해 스케일링 됨
② 그 뒤 데이터는 total cellular read count와 mitochondrial read count의 회귀를 통해 z-score로 스케일링 됨
⑹ ScaleData
① pbmc@assays$RNA@data로부터 pbmc@assays$RNA@scale.data를 생성
② 연산 1. pbmc@assays$RNA@data에 저장돼 있는 normalized expression의 영점을 조정
○ pbmc@assays$RNA@scale.data의 최솟값은 음수가 되는 경우가 대부분
③ 연산 2. 스케일링을 했을 때 값이 10이 초과되는 경우 값을 10으로 할당
○ 이런 할당 연산이 없으면 pbmc@assays$RNA@data와 pbmc@assays$RNA@scale.data 간 spearman correlation 결과는 1 : rank가 보존되므로
④ 진정한 의미의 유전자 발현은 pbmc@assays$RNA@scale.data로 보는 것이 타당함
○ 실제로 TCGA(The Cancer Genome Atlas) 등 bulk-RNA seq에서 쓰는 0-10 scale은 scale.data인 것으로 필자는 추정
⑺ RunPCA
① PCA(principal component analysis)로 차원 축소를 진행 : 클러스터링의 효율을 높이기 위함
② 50개의 주성분이 클러스터링에 활용됨
⑻ JackStraw
① dimensionality를 결정하기 위함으로 이미 레퍼런스 상에 어떤 파라미터를 쓰라고 정해줬으면 그 부분을 생략할 수 있음
② num.replicate : 이 값이 1이어도 추후 진행에 차이가 없음
⑼ ScoreJackStraw
① dimensionality를 결정하기 위함으로 이미 레퍼런스 상에 어떤 파라미터를 쓰라고 정해줬으면 그 부분을 생략할 수 있음
⑽ FindClusters를 통해 클러스터링을 진행하는데 일반적으로 0.5의 resolution이 사용됨
① resolution : 이 값이 커야 cell type을 더 다양하게 분류해줌
② FindClusters는 graph-based approach를 사용함
⑾ RunUMAP
① dims : 맵의 개형을 다소 바꿀 수 있음
⑿ FindAllMarkers에서 Wilcoxon rank sum test가 사용됨
① 각 클러스터 내 25% 이상의 cell에서 발현하고 있고 fold change가 0.25 (log scale)보다 커야 marker gene으로 취급됨
⒀ DotPlot
① https://github.com/satijalab/seurat/issues/2491
② 위 코드에는 안 나와있지만 상당히 자주 사용함
⒁ FeaturePlot
① https://github.com/satijalab/seurat/issues/2400
② 위 코드에는 안 나와있지만 상당히 자주 사용함
2. 커스토마이징(customizing) [목차]
⑴ 개요 : 자기 데이터를 Seurat으로 분석하고 싶은 경우
⑵ filtered_gene_bc_matrices/hg19/ 디렉터리 안에는 총 세 개의 파일이 저장돼 있음
① barcodes.tsv : 각 single cell들의 바코드 네임
② genes.tsv : 인간 유전자의 정식 명칭과 별명
③ matrix.mtx : 각 바코드와 유전자에서 그 유전자의 발현량을 sparse matrix로 나타낸 것
⑶ 1st. matrix.mtx를 수정
① 1st - 1st. 데이터만 표시된 행렬이 주어질 것 (단, rownames와 colnames가 있을 것)
② 1st - 2nd. 0이 아닌 각 데이터가 나타내는 유전자가 genes.tsv에서 몇 번째인지를 첫 번째 값으로 지정
③ 1st - 3rd. 0이 아닌 각 데이터가 나타내는 single cell의 바코드가 무엇인지를 두 번째 값으로 지정
④ 1st - 4th. 0이 아닌 각 데이터에 표시된 발현량을 세 번째 값으로 지정
⑤ 1st - 5th. 첫 번째, 두 번째, 세 번째 값을 한 행으로 함 → 행들을 쌓아서 전체 행렬을 구성 (sparse matrix 구성)
⑥ 1st - 6th. 전체 행렬을 텍스트 파일로 저장하고 .mtx 파일로 변환
⑦ 1st - 7th. 기존 디렉터리에 있던 matrix.mtx 파일과 맞바꾸고 다음과 같이 첫 세 줄을 추가함
Figure. 2. matrix.mtx의 형태
%%MatrixMarket matrix coordinate real general
%
32738 1047 2572919
○ 32,738 : 인간 유전자의 개수로 변경하지 않아도 됨
○ 1,047은 바코드의 개수로 실험에 따라 변경해야 하는 값임
○ 2,572,919는 0이 아닌 데이터의 개수로 실험에 따라 변경해야 하는 값임
⑧ matrix.mtx보다는 matrix.mtx.gz로 read하는 게 권장됨
⑷ 2nd. barcodes.tsv를 수정 : 실험 데이터가 됐던 single cell들의 개수만큼의 바코드를 아무렇게나 부여 (단, 겹치지 않을 것)
① matrix.mtx에 표기된 바코드 순서가 barcodes.tsv의 바코드 개수의 범위를 넘어가면 다음 오류 메시지를 출력함
readMM(): column values 'j' are not in 1:nc
② barcodes.tsv보다는 barcodes.tsv.gz로 read하는 게 권장됨
⑸ 3rd. genes.tsv는 수정하지 않아도 됨
① Seurat 버전이 최신화됨에 따라 genes.tsv를 features.tsv로 바꿔주어야 함
② genes.tsv보다는 genes.tsv.gz로 read하는 게 권장됨
3. 트러블슈팅 [목차]
⑴ "파일이 너무 커서 메모장에서 열 수 없습니다."
① 텍스트 파일이 1 GB 근처를 넘기 시작하면 메모장을 통해 위와 같은 내용을 추가하기 불가능함
② 워드패드는 너무 오래 걸리는 데다 제대로 원하는 기능을 수행하지도 못함
③ gVim 8.0 (윈도우), sublime txt (맥북)등과 같은 프로그램을 통해 용량이 상당히 큰 텍스트 파일을 편집할 수 있음 (ref)
⑵ "Directory provided does not exist"
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
# Error in Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/") :
# Directory provided does not exist
① 파일을 직접 다운받아서 압축을 푼 뒤 생성된 문서의 디렉토리 주소로 바꿔주어야 함
② 파일은 여러 사이트에 공개돼 있음 (ref)
③ 다음과 같이 코드를 수정
pbmc.data <- Read10X(data.dir = "C:/Users/sun/Desktop/filtered_gene_bc_matrices/hg19/")
⑶ Barcode file missing. Expecting barcodes.tsv.gz
① 최근 Seurat 파이프라인은 구버전과 달리 barcodes.tsv, genes.tsv, matrix.mtx가 아니라 barcodes.tsv.gz, features.tsv.gz, matrix.mtx.gz 파일을 요구하고 있음
② 이렇게 구버전 Seurat에 맞춰진 데이터셋을 읽고자 하는 경우 Read10X 대신 ReadMtx를 사용해야 함
⑷ "에러: package or namespace load failed for ‘Seurat’ in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()) versionCheck = vI[[j]]): ‘multtest’이라고 불리는 패키지가 없습니다"
install.packages("BiocManager")
BiocManager::install('multtest')
install.packages("Seurat")
library(Seurat)
⑸ "에러: 크기가 13.3 GB인 벡터를 할당할 수 없습니다"
① 원인 : 프로세스가 메모리 한계를 초과한 경우
○ 일반적으로 pbmc ← ScaleData(object = pbmc, features = all.genes) 부분에서 메모리 수요가 큼
○ pbmc ← JackStraw(object = pbmc, num.replicate = 100) 부분에서도 수요가 큰데 num.replicate = 1로 바꾸면 됨
○ memory.limit으로 늘릴 수 있는 메모리는 PC의 램의 메모리를 초과할 수 없음
○ 작업관리자 → 성능 → 메모리
○ 총 메모리의 용량이 얼마인지 알 수 있음
○ 기본적으로 Seurat를 돌리는 데 7 GB를 잡아먹고 데이터 크기에 따라 추가적인 메모리를 요구함 (e.g., 13.3 GB)
○ 사용된 슬롯을 확인하고 추가로 램을 구입하는 게 하나의 방법이 될 수 있음
○ (주석) 데이터 분석을 한다고 하면 램 32 GB는 구비를 해야 할 듯
② 해결방법 (R 4.1.x 버전 혹은 그 이하)
memory.size(max = TRUE) # OS로부터 R이 사용 가능한 메모리
memory.size(max = FALSE) # 현재 사용중인 메모리 크기
memory.limit(size = NA) # 컴퓨터의 최대 메모리 한계치
memory.limit(size = 50000) # 컴퓨터의 최대 메모리 한계치 약 49GB로 높이기
③ 해결방법 (R 4.2.x 이상 버전)
○ "memory.limit() is no longer supported"라는 에러 메시지가 발생
○ 이 경우 다음과 같이 R의 최대 메모리 한계치를 늘릴 수 있음 (레퍼런스)
library(usethis)
usethis::edit_r_environ()
# Then, type "R_MAX_VSIZE=100Gb" and save .Renviron file.
○ 다만, R 4.2.x 이상 버전에서는 자동으로 메모리를 늘리기 때문에 별도로 걱정할 필요는 없음 (레퍼런스)
⑹ "dimnamesGets(x, value)에서 다음과 같은 에러가 발생했습니다: invalid dimnames given for “dgTMatrix” object"
① 이 문제는 위와 같이 genes.tsv에서 값들이 잘 정렬되지 않아 Read10X 부분에서 문제를 일으키는 것
② 모든 tuple들이 3열이 되도록 잘 정리한다면 다음 단계가 smooth하게 진행될 것임
⑺ An old seurat object 23341 genes across 624 samples
① 문제점 : old seurat object는 내용을 들여다 보거나 분석을 진행하기 어려움
② 해결방법 : UpdateSeuratObject를 사용
load('~/Downloads/droplet_Heart_and_Aorta_seurat_tiss.Robj')
pbmc = UpdateSeuratObject(tiss)
DimPlot(pbmc)
⑻ Seurat object로부터 matrix.mtx, barcodes.tsv, genes.tsv 파일을 저장하는 방법 (R)
library(Seurat)
counts <- seurat_obj@assays$RNA@counts
barcodes <- rownames(seurat_obj@meta.data)
features <- rownames(counts)
library(Matrix)
writeMM(counts, file = "~/Downloads/matrix.mtx")
write.table(barcodes, file = "~/Downloads/barcodes.tsv", row.names = FALSE, col.names = FALSE, quote = FALSE, sep = "\t")
gene_ids <- features
gene_names <- features
features_df <- data.frame(gene_ids, gene_names)
write.table(features_df, file = "~/Downloads/features.tsv", row.names = FALSE, col.names = FALSE, quote = FALSE, sep = "\t")
⑼ matrix.mtx, barcodes.tsv, features.tsv (or genes.tsv)가 있을 때 h5ad 생성하기 (Python) : 단, 세 파일 모두 prefix가 동일하다고 가정
import os
import glob
from scipy.io import mmread
import pandas as pd
import scanpy as sc
# Set the target directory here
target_directory = '/path/to/your/directory'
# Iterate over all .mtx files in the target directory
for mtx_file in glob.glob(os.path.join(target_directory, '*matrix.mtx')):
prefix = mtx_file[:-10] # Remove 'matrix.mtx' to get the prefix
# Construct file names for barcodes and genes
barcodes_file = prefix + 'barcodes.tsv'
features_file = prefix + 'features.tsv'
genes_file = prefix + 'features.tsv' if os.path.exists(prefix + 'features.tsv') else prefix + 'genes.tsv'
# Check if all files exist
if os.path.exists(mtx_file) and os.path.exists(barcodes_file) and os.path.exists(genes_file):
# Load the data
matrix = mmread(mtx_file).T.tocsr() # Transpose to have genes as columns
genes = pd.read_csv(genes_file, header=None, sep='\t')
barcodes = pd.read_csv(barcodes_file, header=None, sep='\t')
# Create an AnnData object
adata = sc.AnnData(X=matrix, obs={'barcodes': barcodes[0].values}, var={'gene_ids': genes[0].values, 'gene_symbols': genes[1].values})
# Output file name
h5ad_file = prefix + '.h5ad'
# Save the AnnData object
adata.write(h5ad_file)
# Delete the original files
os.remove(mtx_file)
os.remove(barcodes_file)
os.remove(genes_file)
print("Conversion complete. Original files have been deleted.")
⑽ tissue_dir 디렉토리 내 matrix.mtx, barcodes.tsv, features.tsv 및 spatial 폴더가 있을 때, Visium 데이터를 Python 상에서 읽는 코드
import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.image import imread
import numpy as np
import cv2
import math
from pathlib import Path, PurePath
from typing import Union, Dict, Optional, Tuple, BinaryIO
import json
def to_spatial(adata, path, library_id = 'TNBC1'):
adata.uns["spatial"][library_id] = dict()
path = Path(path)
files = dict(
tissue_positions_file=path / 'spatial/tissue_positions_list.csv',
scalefactors_json_file=path / 'spatial/scalefactors_json.json',
hires_image=path / 'spatial/tissue_hires_image.png',
lowres_image=path / 'spatial/tissue_lowres_image.png',
)
# check if files exists, continue if images are missing
for f in files.values():
if not f.exists():
if any(x in str(f) for x in ["hires_image", "lowres_image"]):
logg.warning(
f"You seem to be missing an image file.\n"
f"Could not find '{f}'."
)
else:
raise OSError(f"Could not find '{f}'")
adata.uns["spatial"][library_id]['images'] = dict()
for res in ['hires', 'lowres']:
try:
adata.uns["spatial"][library_id]['images'][res] = imread(
str(files[f'{res}_image'])
)
except Exception:
raise OSError(f"Could not find '{res}_image'")
# read json scalefactors
adata.uns["spatial"][library_id]['scalefactors'] = json.loads(
files['scalefactors_json_file'].read_bytes()
)
#adata.uns["spatial"][library_id]["metadata"] = {
# k: (str(attrs[k], "utf-8") if isinstance(attrs[k], bytes) else attrs[k])
# for k in ("chemistry_description", "software_version")
# if k in attrs
#}
# read coordinates
positions = pd.read_csv(files['tissue_positions_file'], header=None)
positions.columns = [
'barcode',
'in_tissue',
'array_row',
'array_col',
'pxl_col_in_fullres',
'pxl_row_in_fullres',
]
positions.index = positions['barcode']
adata.obs = adata.obs.join(positions, how="left")
adata.obsm['spatial'] = adata.obs[
['pxl_row_in_fullres', 'pxl_col_in_fullres']
].to_numpy()
adata.obs.drop(
columns=['barcode', 'pxl_row_in_fullres', 'pxl_col_in_fullres'],
inplace=True,
)
return adata
def load_adata(tissue_dir):
adata1 = sc.read_mtx(tissue_dir + 'matrix.mtx')
adata1 = adata1.transpose()
cellname = pd.read_csv(tissue_dir + 'barcodes.tsv', sep="\t", header=None, index_col=1, names=['idx', 'barcode'])
cellname.index = cellname['idx']
featurename = pd.read_csv(tissue_dir + 'features.tsv', sep='\t', header=None, index_col=1, names=['gene_ids', 'feature_types'])
adata1.var = featurename
adata1.obs = cellname
adata1.uns["spatial"] = dict()
adata1 = to_spatial(adata1, tissue_dir, library_id='Library_1')
adata1 = adata1[adata1.obs.in_tissue == 1,:]
adata1.var_names_make_unique()
adata1.obs['unnormalized_counts'] = adata1.X.sum(axis=1).A1
sc.pp.normalize_total(adata1, inplace=True)
adata1.obs['normalized_counts'] = adata1.X.sum(axis=1).A1
sc.pp.log1p(adata1)
sc.pp.highly_variable_genes(adata1, flavor="seurat", n_top_genes=2000)
return adata1
⑾ matrix.mtx.gz, barcodes.tsv, 또는 features.tsv가 깨진 파일인 경우
① matrix.mtx.gz 파일이 깨진 경우 R 환경에서 readMM으로 읽어들인 뒤 writeMM으로 깨진 파일을 복구
② barcodes.tsv 및 features.tsv의 경우 R 환경에서 read.table로 읽어들인 뒤 write.table로 깨진 파일을 복구
⑿ 똑같은 코드를 쓰면 똑같은 결과가 나오지만 컴퓨터 종류에 따라 결과가 다르게 나올 때도 있고 같게 나올 때도 있음
① 추정하는 이유 : 초기 조건, 즉 seed가 약간 차이가 나는 것으로 보임
⒀ 입력 파일이 .rds 파일인 경우
# when loading .rds file
pbmc <- readRDS("Example.rds")
# when saving .rds file
saveRDS(pbmc, file = "Example.rds")
## You can even plot pbmc, when the .rds file is fully processed file
DimPlot(object = pbmc, reduction = "umap")
① .rds 파일 : R에서 하나의 변수를 저장하고 있음
⒁ 입력 파일이 .RData 파일 또는 .Rdata 파일인 경우
# when loading .RData file
load(file="Example.RData")
# when saving .RData file
## method 1.
save.image(file="Example.RData")
## method 2.
save(x, file='Example.Rdata')
① .RData 파일 : R에서 모든 변수를 저장하고 있음
⒂ 입력 파일이 .RDS 파일 또는 .Rds 파일인 경우
# when loading .RDS file
## method 1.
load(file="Example.RDS")
## method 2.
pbmc <- readRDS('Example.RDS')
① .RDS 파일 : R에서 모든 변수 혹은 단일 변수를 저장하고 있음
⒃ 입력 파일이 .h5 파일인 경우
① R : 다음 코드는 .h5 파일을 sparse matrix로 읽어들임
library(Seurat)
data = Read10X_h5("data.h5")
② 파이썬 : 다음 코드는 .h5 파일을 Anndata object로 읽어들임
import scanpy as sc
data = sc.read_10x_h5('data.h5')
⒄ 입력 파일이 .h5ad 파일인 경우
① R : 다음 코드는 .h5ad 파일을 Anndata object로 읽어들임 (2023년 3월 신설)
library(anndata)
data = read_h5ad("data.h5ad")
② 파이썬 : 다음 코드는 .h5 파일을 Anndata object로 읽어들임
pbmc = scanpy.read_h5ad("data.h5ad")
③ Seurat 객체 pbmc3k를 h5ad로 저장하기 (레퍼런스)
library(Seurat)
library(SeuratData)
library(SeuratDisk)
SaveH5Seurat(pbmc3k, filename = "~/Downloads/pbmc3k.h5Seurat")
Convert("~/Downloads/pbmc3k.h5Seurat", dest = "h5ad")
⒅ 입력 파일이 .rda 파일인 경우
load("/Users/parkjeongbin/Downloads/ifnb.SeuratData/data/ifnb.rda")
① .rda 파일 : R에서 하나의 변수를 저장하고 있음. .rds 파일과 달리 return 타입이 없는 load 함수를 사용하여 변수를 로드함
② 참고로 위 ipynb.rda 파일은 다음과 같이 RStudio 상에서 설치하였음
library(SeuratData)
InstallData("ifnb")
⒆ 입력 파일이 .Rmd 파일인 경우
① .ipynb 파일과 .Rmd 파일은 상호 변환 가능
② .Rmd 파일에서 .ipynb 파일을 생성하는 코드 (레퍼런스)
devtools::install_github("mkearney/rmd2jupyter")
library(rmd2jupyter)
rmd2jupyter("MyFILE.Rmd")
⒇ 입력 파일이 .Robj 파일인 경우 (ref)
① load() 함수로 입력 가능
load('Data.Robj')
⒇ 입력 파일이 .cel 파일인 경우 (ref1, ref2)
① .cel file은 Affymetrix로 생성된 raw data 파일 확장자
② 코드
b <- ReadAffy(filenames='/Users/parkjeongbin/Desktop/614615111406.E02.CEL')
#Extract the perfect match probe-level intensities
dim(pm(b)) #[1] 247965 1
#Phenotypic data : sample information 6X17
dim(pData(b)) #[1] 1 1
#Plateform used for gene information
annotation(b) #[1] "hgu133a"
#Retrieving other types of information
## Use the cfdName() method on the AffyBatch
cdfName(b) #[1] "HG-U133A"
## Use the featureNames() method on the AffyBatch
featureNames(b)
## Use the length() method to count the number of items in a vector
length(featureNames(b)) #[1] 22283
## Use the length() method to count the number of items in the vector containing names of all probes
length(probeNames(b)) #[1] 247965
#Normalize
data.rma = rma(b)
data.matrix = exprs(data.rma)
dim(data.matrix) #[1] 22283 1
#Convert Affymetrix probe ID to ensembl_gene_id, gene_name
## https://www.biostars.org/p/328065/#328328
## https://www.biostars.org/p/332461/#332474
BiocManager::install("biomaRt")
library(biomaRt)
dat<-rownames(data.matrix)
require("biomaRt")
mart <- useMart("ENSEMBL_MART_ENSEMBL")
mart <- useDataset("hsapiens_gene_ensembl", mart)
annotLookup <- getBM(
mart = mart,
attributes = c(
"affy_hg_u133_plus_2",
# "ensembl_gene_id",
# "gene_biotype",
"external_gene_name"),
filter = "affy_hg_u133_plus_2",
values = dat,
uniqueRows=TRUE)
○ 위 예제의 특징은 다음과 같음
○ single file을 읽어들여 sample이 1개밖에 없음 : 그러나 폴더 자체를 읽어들일 수도 있음
○ HG-U133A라는 assay를 이용함
○ 22283개의 feature가 있음 : feature는 곧 gene, probe_id를 의미함
○ gene expression은 최대 14로 scaling 됨
⒇ 소스 파일이 .ipynb 파일인 경우
① jupyter notebook 파일로 코드와 코드가 실행되는 과정이 기록돼 있음
② Python 뿐만 아니라 R로도 .ipynb 파일을 생성할 수 있음
③ 내용을 확인하는 방법 : Anaconda, Google Colab, GitHub
⒇ 소스 파일이 .Rmd 파일인 경우
① Rmd 파일을 ipynb 파일로 변환시키는 패키지
### conver Rmd into ipynb ###
devtools::install_github("mkearney/rmd2jupyter")
library(rmd2jupyter)
rmd2jupyter("filename.Rmd")
# Rmd 파일이 있는 path에 finename.ipynb 생성
입력: 2019.11.25 20:22
수정: 2022.12.28 23:39
'▶ 자연과학 > ▷ 생물정보학' 카테고리의 다른 글
| 【생물정보학】 MIA 분석의 이해 및 실행 (0) | 2021.01.02 |
|---|---|
| 【생물정보학】 scater로 cell type 결정하기 (0) | 2019.12.20 |
| 【생물정보학】 Cell Type Classification Pipeline (0) | 2019.11.22 |
| 【생물정보학】 TCGA DATA 얻는 법 (2) | 2019.08.26 |
| 【생물정보학】 생물정보학 노트 (6) | 2018.02.10 |
최근댓글