هر سلول از بدن ما ۲۰ هزار ژن یکسان دارد، و این ژن‌ها در سلول‌های مختلف و در شرایط مختلف ممکن است روشن یا خاموش باشند (transcription). در واقع عاملی که باعث تفاوت سلول‌های پوست، با سلول‌های عصب می‌شود، همین تفاوت الگوی فعال بودن دسته‌های مختلف ژن‌ها، در این سلول‌هاست.
با استفاده از تکنولوژی سنجش بیان ژن‌ها در تک‌سلول‌ها (single-cell RNA-Seq) می‌توان نسبت به پیچیدگی سامانه‌های زیستی متشکل از انواع سلول‌ها دید پیدا کرد. به عنوان نمونه از این تکنولوژی برای آگاهی پیدا کردن به گوناگونی سلول‌های سرطان در تومور و فرآیند‌های دارای پتانسیل برای مداخله دارویی به منظور مهار تومور دید پیدا کرد. به عنوان نمونه‌ای دیگر، از این تکنولوژی برای مطالعه کارکرد انواع سلول‌های ایمنی و همچنین مهندسی سلول‌های ایمنی به منظور مهار تومور استفاده می‌شود.
در این کارگاه در پژوهشگاه رویان یک workflow استاندارد برای آنالیز داده‌های single-cell RNA-Seq یک دیتاست عموما منتشر شده، ارائه کردم. طی چند ماه بعد از این کارگاه ایمیل‌هایی مبنی بر به اشتراک گذاری مطالب این کارگاه دریافت کردم، و تصمیم گرفتم این مطلب آموزشی را در اینجا منتشر کنم. در این مطلب سعی میکنم اول-تا-آخر مراحل آنالیز داده‌های منتشر شده در این مقاله را مطرح کنم.
https://elifesciences.org/articles/55320
در مراحل مختلف این مطلب، نحوه رسیدن به شکل‌های منتشر در این مقاله روشن می‌شود.
این مطالب آموزشی در رابطه به تکنولوژی 10x Genomicsنوشته شده است، که یکی از پلتفورم‌های پرطرفدار برای کاربر single-cell RNA-Seq است. خروجی این پلتفورم libraryهای آماده ورود به ماشین توالی‌یابی Illumina است. این مطلب آموزشی، مانند بسیاری از ابزار بیوانفورماتیک، با فرض محیط سیستم عامل linux تالیف شده است.
داده‌های ورودی
نقطه شروع چنین آنالیزی، فایل‌های fastq مربوط به سلول‌ها/شرایط آزمایش مختلف هستند. این فایل‌ها قطعه‌های کوتاه توالی cDNA که مشخص کننده‌ی یک کپی بیان شده از RNA مربوط به یک ژن مشخص در یک سلول مشخص است را در بر دارند. ساختار این توالی‌های کوتاه در شکل زیر آمده است.
در شکل بالا قسمت قرمز رنگ یک توالی تصادفی است، که برای همه مولکول‌های cDNA حاصل از یک سلول، مشترک است. پس با توجه به این قسمت، توالی‌های مربوط به هزاران سلول حاضر در هر نمونه آزمایش را می‌توان به سلول مربوطه متناظر کرد. همچنین قسمت سبز رنگ مشخص می‌کند این توالی مربوطه به کدام نمونه (نوع سلول/ شرایط آزمایش/ تکرار) است. لازم به یادآوری است که همه توالی‌های مربوط به سلول‌های مختلف و نمونه‌های مختلف ممکن است به همراه هم وارد ماشین Illumina بشوند، و از این تگ‌های توالی می‌توان نمونه‌ها و سلول‌ها را از هم تفکیک کرد. همچنین با بررسی توالی قسمت insert و مقایسه آن با ژنوم مرجع (مثلا hg38) می‌توان تعیین کرد این cDNA مربوط به کدام ژن است.
نحوه دریافت داده‌های خام
نقطه شروع ما فایل‌های fastq هستند. برای یک آزمایش خاص که هدف این مطلب است این فایل‌ها در لینک زیر در دامنه عمومی قرار گرفته است، و مربوط به دیتاستی است که به همراه مقاله مربوطه منتشر شده است.
https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA605000&o=acc_s%3Aa
اما این اینترفیس حاضر در لینک فقط اجازه دریافت فایل از طریق سرویس‌های ابری AWS یا GCP را می‌دهد. برای دانلود فایل‌ها به صورت محلی، می‌توان از ابزار sra-toolkit می‌توان این داده را به صورت محلی دانلود کرد. این ابزار به عنوان آرگومان ورودی شناسه SRA دیتاست هدف است و خروجی آن فایل‌های fastq است که دانلود می‌شود.
مراحل نصب sra-toolkit اینجا شرح داده شده است.
نحوه استفاده از sra-toolkit برای دریافت داده‌های مربوط به یک دیتاست مشخص با شناسه SRA اینجا شرح داده شده است. مطابق این دستورالعمل با دستور‌های زیر می‌توان فایل‌های fastq این دیتاست مشخص را دریافت کرد. توجه داشته باشید که مجموع حجم همه نمونه‌های موجود در این دیتا ست ۱۴۰ گیگابایت است. پس از نصب این ابزار، فایل قابل اجرا با نام fasterq-dump ایجاد می‌شود و به صورت زیر از آن استفاده می‌کنیم.
$ fasterq-dump --split-files PRJNA605000
نصب ابزار مورد نیاز
از cellranger برای شمارش تعداد ترانسکریپت مربوط به هر ژن در هر سلول بر اساس داده‌های fastq استفاده می‌کنیم. این ابزار در محیط command-line لینوکس کار می‌کند. از لینک زیر می‌توانید درخواست لینک دانلود این ابزار را بدهید.
https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest
با استفاده از دستور زیر می‌توانید cellranger را از لینک تولید شده بالا را دانلود کنید.
curl -o cellranger-4.0.0.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-4.0.0.tar.gz?Expires=BLA..BLA!!GET_YOUR_OWN_LINK!!!!
این فایل 1GB حجم دارد. اگر curl در محیط شما نصب نیست، می‌توانید با دستور زیر آن‌را نصب کنید
$ sudo apt-get install curl
با دستور زیر فایل زیپ شده cellranger را آنزیپ و نصب می‌کنیم.
$ tar -xvzf cellranger-4.0.0.tar.gz
نرم‌افزار در دیرکتوری با نام cellranger-4.0.0 قرار می‌گیرد.
ادامه مراحل را به زبان انگلیسی می‌آورم چون از قبل به این صورت مستند کرده بودم. امیدوارم قابل استفاده باشد، و بتوانم سر فرصت آن را به صورت مناسبی به فارسی برگردانم.
################################
Running cellranger mkref
################################
#Next these fastq files need to mapped to the reference genome. To do this, we primarily need to prepare the reference genome in the particular format needed for cellranger. This is done using the 'cellranger mkref' command. The help of this command can be consulted with th following command:
/home/abbas/myProjects/scRNAseq_data_analysis/run_cellranger_count/glu12h$ ~/myProjects/scRNAseq_data_analysis/cellranger-4.0.0/cellranger mkref --help
#In order to run the mkref command to input files are needed: I) a fasta file that contains the genomic DNA sequence for each chromosome. II) A gtf file that annotates the feature (e.g. genes) loci on the reference genome. Both of these files were downloaded from UCSC genome browser. The fasta file was downloaded from the 'Sequence and Annotation Downloads' section of this website, and the gtf file was downloaded from the 'Table Browser' section.
#We can then execute the mkref command like below, assuming you have the fasta and gtf files in the current directory:
/home/abbas/myProjects/scRNAseq_data_analysis/cellranger-4.0.0/cellranger mkref --genome=S288C --fasta=sacCer3.fa --genes=SacCer3.gtf
################################
Running cellranger count
################################
#Afterwards, we proceed to run the 'cellranger count' command. This command maps the reads to the genomes, and creates the bam files, and as well carries out the UMI counting for each gene. We need to run this command for each of the samples separately, in separate folders. You can run this commands in parallel in separate terminal windows. By inspecting the help of this command, we can see the major inputs to this command are, the fastq files corresponding to this sample, from both NextSeq and HiSeq runs, as well as the reference genome that we just created. Moreover some metadata, such as the expected cell counts, and the chemistry version of the kit used are provided. Finally, you can put limits on memory and cpu usage of this command. This is the most computationally intensive step of the analysis.
#For the glu12h sample, this command was executed like this:
mkdir glu12h
cd glu12h
/home/abbas/myProjects/scRNAseq_data_analysis/cellranger-4.0.0/cellranger count --id glu12h --fastqs /media/abbas/Abbas_5TB/scRNAseq_data_analysis/run_cellranger_mkfastq/mk_fastq_HiSeq/HiSeqMkFastqOut/outs/fastq_path/HJ5FTBBXX/glu12h,/media/abbas/Abbas_5TB/scRNAseq_data_analysis/run_cellranger_mkfastq/mk_fastq_NextSeq/NextSeqMkFastqOut/outs/fastq_path/H27JJBGX3/glu12h --sample glu12h --transcriptome /home/abbas/myProjects/scRNAseq_data_analysis/run_cellranger_mkref/from_UCSC/S288C --expect-cells 2070 --localcores 15 --localmem 40 --chemistry SC3Pv2
#Similar command was executed for the four remaining samples (glu6h, inlag1h, inlag3h, mix_glu_mal) in parallel. On a laptop with a Core i7-8565U CPU, and 32GB of RAM, these commands took about 30 hours to run.
################################
Running cellranger aggr
################################
The final step in the cellranger pipeline is the 'cellranger aggr' command, which aggregates the counts from the 5 samples. In order to run this command we need an aggregation_input.csv file which has one row per sample, and includes the path to the output of the 'cellranger count' command. This the content of such csv file:
library_id,molecule_h5
glu12h,/media/abbas/scRNAseq_data_analysis/run_cellranger_count/glu12h/glu12h/outs/molecule_info.h5
glu6h,/media/abbas/scRNAseq_data_analysis/run_cellranger_count/glu6h/glu6h/outs/molecule_info.h5
inlag1h,/media/abbas/scRNAseq_data_analysis/run_cellranger_count/inlag1h/inlag1h/outs/molecule_info.h5
inlag3h,/media/abbas/scRNAseq_data_analysis/run_cellranger_count/inlag3h/inlag3h/outs/molecule_info.h5
mix_glu_mal,/media/abbas/scRNAseq_data_analysis/run_cellranger_count/mix_glu_mal/mix_glu_mal/outs/molecule_info.h5
#We the run the aggr command like this
/home/abbas/myProjects/scRNAseq_data_analysis/cellranger-4.0.0/cellranger aggr --id=AGG_all_ncellsAuto --csv=aggregation_input.csv
#The output of this step can be used as input to Seurat R package for secondary analysis
##################################
################################
Secondary analysis in Seurat
################################
#To use Seurat you must primarily have R and Rstudio installed. (https://www.r-project.org/, https://rstudio.com/)
#You then need to install the Seurat package. To do this run the following command in the R enivornment:
install.packages('Seurat')
#ggplot2 and dplyr are two additional R packages that we use here:
install.packages('ggplot2')
install.packages('dplyr')
#We then use the Read10X function of the Seurat pacakge, in order to read in the output generated in the last step of cellranger analysis:
data_merge <-
Read10X(data.dir = "/media/abbas/scRNAseq_data_analysis/run_cellranger_aggr/AGG_all_ncellsAuto/outs/filtered_feature_bc_matrix/")
#Inspect the data matrix that you just read in, using the commands below, can you understand the logic behind column names and row names?
print(colnames(data_merge))
print(rownames(data_merge))
print(dim(data_merge))
#Each column in this matrix corresponds to a single cell, and is denoted by a molecular barcode, followed by a numerical index (-1,-2,...,-5). This numerical index describes which of the 5 samples, each particular cell comes from.
#Each row in this matrix corresponds to a gene.
#We then create a Seurat object from this matrix:
data_merge_obj <- CreateSeuratObject(counts = data_merge)
#assigning sample names, based on the numerical index
data_merge_obj@meta.data$timepoint <- 'na'
data_merge_obj@meta.data[grepl('-1$',colnames(data_merge)),'timepoint'] <- 'glu12h'
data_merge_obj@meta.data[grepl('-2$',colnames(data_merge)),'timepoint'] <- 'glu6h'
data_merge_obj@meta.data[grepl('-3$',colnames(data_merge)),'timepoint'] <- 'lag1h'
data_merge_obj@meta.data[grepl('-4$',colnames(data_merge)),'timepoint'] <- 'lag3h'
data_merge_obj@meta.data[grepl('-5$',colnames(data_merge)),'timepoint'] <- 'mixglumal'
#In S cerevisiae, mitochondrial gene names start with Q.
#We need this, since the percentage of mitochondrial genes is a common QC metric.
data_merge_obj[["percent.mt"]] <- PercentageFeatureSet(data_merge_obj, pattern = "^Q")
#plot three major qc metrics
p1 <- VlnPlot(data_merge_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
ggsave('qc_metrics.png', p1)
#apply the filtering
data_merge_obj <-
data_merge_obj[,(grepl('-1$', colnames(data_merge_obj)) & data_merge_obj$nFeature_RNA < 2000 & data_merge_obj$percent.mt < 0.2) |
(grepl('-2$', colnames(data_merge_obj)) & data_merge_obj$nFeature_RNA < 2000 & data_merge_obj$percent.mt < 0.2) |
(grepl('-5$', colnames(data_merge_obj)) & data_merge_obj$nFeature_RNA < 2000 & data_merge_obj$percent.mt < 0.2) |
(grepl('-3$', colnames(data_merge_obj)) & data_merge_obj$nFeature_RNA < 1500 & data_merge_obj$percent.mt < 0.5) |
(grepl('-4$', colnames(data_merge_obj)) & data_merge_obj$nFeature_RNA < 1500 & data_merge_obj$percent.mt < 1.5) ]
#############################
#checking for possible batch effects as described in the paper, using glu12h and mix_glu_mal samples. The key point here is that half of the cells in the mix_glu_mal sample come from the same condition as glu12h.
glu12h_mixglumal_obj <- data_merge_obj[,grepl('-5$|-1$', colnames(data_merge_obj))]
glu12h_mixglumal_obj <- NormalizeData(glu12h_mixglumal_obj)#, normalization.method = 'CLR', margin = 1)
glu12h_mixglumal_obj <- FindVariableFeatures(glu12h_mixglumal_obj, selection.method = "vst", nfeatures = 2000)
glu12h_mixglumal_obj <- ScaleData(glu12h_mixglumal_obj, verbose = FALSE)
glu12h_mixglumal_obj <- RunPCA(glu12h_mixglumal_obj, npcs = 20, verbose = FALSE)
glu12h_mixglumal_obj <- RunUMAP(glu12h_mixglumal_obj, reduction = "pca", dims = 1:20 , n.epochs = 1000)
#color pallet
cbp1 <- c("#E69F00", "#56B4E9", "#009E73",
"#F0E442", "#0072B2", "#D55E00", "#CC79A7")
p2 <-
DimPlot(glu12h_mixglumal_obj, reduction = "umap", group.by = "timepoint")+
scale_colour_manual(values=cbp1) +
theme(legend.position = 'none') +
xlab('') + ylab('')
ggsave('batchEffect_glu12h_mixglumal_1911214.png', p2, dpi=250, width = 5, height = 5)
###############################################################
#plotting the expression of maltose growth markers in the cells
p_ls <-
FeaturePlot( glu12h_mixglumal_obj, c('MAL11','MAL12','MAL31', 'MAL32', 'IMA1', 'HXK1'),
order = F, pt.size = 0.5,
max.cutoff = 0.0001, cols = c('gray50','red'),
combine = FALSE)
for(i in 1:length(p_ls)) {
p_ls[[i]] <- p_ls[[i]] + NoLegend() + NoAxes()
}
png('glu12hmixglumal_markers_mal_191214.png', res=600, units = 'in', width = 8, height = 4)
cowplot::plot_grid(plotlist = p_ls)
dev.off()
###############################################
#####################################
#making a umap plot with all samples, separating the maltose positive cells from the mix sampple
mixglumal_obj <- data_merge_obj[,grepl('-5$', colnames(data_merge_obj))]
mixglumal_exprs <- as.matrix(GetAssayData(mixglumal_obj, slot = "counts"))
mixglumal_exprs_mal_markers <- mixglumal_exprs[c('MAL31','MAL32','HXK1','IMA1'),]
cells_MAL <- names(colSums(mixglumal_exprs_mal_markers)[colSums(mixglumal_exprs_mal_markers)>3])
cells_mix_not_MAL <- setdiff(colnames(mixglumal_obj), cells_MAL)
data_wMAL_obj <-
data_merge_obj[,!colnames(data_merge_obj)%in%cells_mix_not_MAL]
data_wMAL_obj@meta.data[grepl('-5$',colnames(data_wMAL_obj)),'timepoint'] <- 'mal'
data_wMAL_obj <- NormalizeData(data_wMAL_obj)#, normalization.method = 'CLR')
data_wMAL_obj <- FindVariableFeatures(data_wMAL_obj, selection.method = "vst", nfeatures = 2000)
data_wMAL_obj <- ScaleData(data_wMAL_obj, verbose = FALSE)
data_wMAL_obj <- RunPCA(data_wMAL_obj, npcs = 20, verbose = FALSE)
data_wMAL_obj <- RunUMAP(data_wMAL_obj, reduction = "pca", dims = 1:20 , n.epochs = 1000)
cbp1 <- c("#E69F00","#CC79A7", "#009E73","gray30", "#0072B2")
p <-
DimPlot(data_wMAL_obj, reduction = "umap", group.by = "timepoint", pt.size = 0.1)+
scale_colour_manual(values=cbp1)
ggsave('all_samples_malSep_191214_wLegend.png',p,dpi=600, width = 5, height = 5)
#############################################
#plotting only the samples during maltose shift
data_merge <-
Read10X(data.dir = "/media/abbas/scRNAseq_data_analysis/run_cellranger_aggr/AGG_all_ncellsAuto/outs/filtered_feature_bc_matrix/")
data_merge_obj <- CreateSeuratObject(counts = data_merge)
data_merge_obj@meta.data$timepoint <- 'na'
data_merge_obj@meta.data[grepl('-1$',colnames(data_merge)),'timepoint'] <- 'glu12h'
data_merge_obj@meta.data[grepl('-2$',colnames(data_merge)),'timepoint'] <- 'glu6h'
data_merge_obj@meta.data[grepl('-3$',colnames(data_merge)),'timepoint'] <- 'lag1h'
data_merge_obj@meta.data[grepl('-4$',colnames(data_merge)),'timepoint'] <- 'lag3h'
data_merge_obj@meta.data[grepl('-5$',colnames(data_merge)),'timepoint'] <- 'mixglumal'
data_merge_obj[["percent.mt"]] <- PercentageFeatureSet(data_merge_obj, pattern = "^Q")
#VlnPlot(data_merge_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
data_merge_obj <-
data_merge_obj[,(grepl('-1$', colnames(data_merge_obj)) & data_merge_obj$nFeature_RNA < 2000 & data_merge_obj$percent.mt < 0.2) |
(grepl('-2$', colnames(data_merge_obj)) & data_merge_obj$nFeature_RNA < 2000 & data_merge_obj$percent.mt < 0.2) |
(grepl('-5$', colnames(data_merge_obj)) & data_merge_obj$nFeature_RNA < 2000 & data_merge_obj$percent.mt < 0.2) |
(grepl('-3$', colnames(data_merge_obj)) & data_merge_obj$nFeature_RNA < 1500 & data_merge_obj$percent.mt < 0.5) |
(grepl('-4$', colnames(data_merge_obj)) & data_merge_obj$nFeature_RNA < 1500 & data_merge_obj$percent.mt < 1.5) ]
data_merge_obj <- data_merge_obj[,grepl('-1$|-3$|-4$', colnames(data_merge_obj))]
data_merge_obj <- NormalizeData(data_merge_obj)#, normalization.method = 'CLR', margin = 1)
data_merge_obj <- FindVariableFeatures(data_merge_obj, selection.method = "vst", nfeatures = 2000)
data_merge_obj <- ScaleData(data_merge_obj, verbose = FALSE)
data_merge_obj <- RunPCA(data_merge_obj, npcs = 20, verbose = FALSE)
data_merge_obj <- RunUMAP(data_merge_obj, reduction = "pca", dims = 1:20 , n.epochs = 1000)#,min.dist = 0.5)
cbp1 <- c("#E69F00", "#009E73","gray30",
"#F0E442", "#0072B2", "#D55E00", "#CC79A7")
p <-
DimPlot(data_merge_obj, reduction = "umap", group.by = "timepoint", pt.size = 0.1)+
scale_colour_manual(values=cbp1) +
#theme(legend.position = 'none') +
xlab('') + ylab('')
ggsave('shift_samples_1911214_wLegend.png', p,dpi=600, width = 5, height = 5)
#########################
#clustering on shift samples
data_merge_obj <- FindNeighbors(data_merge_obj, dims = 1:10)
#The resolution paramter in the command below, determines the number of clusters. We tweaked it to have 4 clusters.
data_merge_obj <- FindClusters(data_merge_obj, resolution = 0.12)
cbp1 <- c("#E69F00", "#009E73","gray30",
"#F0E442", "#0072B2", "#D55E00", "#CC79A7")#
p <-
DimPlot(data_merge_obj, label = FALSE, pt.size = 0.1)+
scale_colour_manual(values=cbp1) +
#theme(legend.position = 'none') +
xlab('') + ylab('')
ggsave('shift_samples_191214_4clusters.png',p, dpi=600, width = 5, height = 5)
########################
#differential expression
diff_exp_res <- FindMarkers(data_merge_obj, ident.1 = 2, ident.2 = 3, logfc.threshold = 0.01)
diff_exp_res <- diff_exp_res[abs(diff_exp_res$avg_logFC)>0.25 & diff_exp_res$p_val_adj<0.05, ]
#write the results to a file
write.csv(diff_exp_res, 'diff_exp_res.csv', row.names = F)
#upload the genes in this file to a Gene Ontology webtool, eg
#http://cbl-gorilla.cs.technion.ac.il/
