admin 管理员组文章数量: 1184232
2023年12月23日发(作者：excel表vlookup公式怎么使用)
("DESeq2")("DESeq2/pairwise")# Function to run DESeq2 and get results for all clusters## x is index of cluster in clusters vector on which to run function## A is the sample group to compare## B is the sample group to compare against (base level)get_dds_resultsAvsB <- function(x, A, B){ cluster_metadata <- metadata[which(metadata$cluster_id == clusters[x]), ] rownames(cluster_metadata) <- cluster_metadata$sample_id counts <- pb[[clusters[x]]] cluster_counts <- (counts[, which(colnames(counts) %in% rownames(cluster_metadata))]) #all(rownames(cluster_metadata) == colnames(cluster_counts)) dds <- DESeqDataSetFromMatrix(cluster_counts, colData = cluster_metadata, design = ~ group_id) # Transform counts for data visualization rld <- rlog(dds, blind=TRUE) # Plot PCA DESeq2::plotPCA(rld,
intgroup = "group_id") ggsave(paste0("results/", clusters[x], "_specific_"))
# Extract the rlog matrix from the object and compute pairwise correlation values
rld_mat <- assay(rld) rld_cor <- cor(rld_mat) # Plot heatmap png(paste0("results/", clusters[x], "_specific_")) pheatmap(rld_cor, annotation = cluster_metadata[, c("group_id"), drop=F]) () # Run DESeq2 differential expression analysis dds <- DESeq(dds) # Plot dispersion estimates png(paste0("results/", clusters[x], "_dispersion_")) plotDispEsts(dds) () # Output results of Wald test for contrast for A vs B contrast <- c("group_id", levels(cluster_metadata$group_id)[A], levels(cluster_metadata$group_id)[B]) # resultsNames(dds) res <- results(dds, contrast = contrast,
alpha = 0.05) res <- lfcShrink(dds, contrast = contrast,
res=res) # Set thresholds padj_cutoff <- 0.05 # Turn the results object into a tibble for use with tidyverse functions res_tbl <- res %>% ()
%>% rownames_to_column(var="gene") %>% as_tibble() (res_tbl, paste0("DESeq2/pairwise/", clusters[x], "_", levels(cluster_metadata$group_id)[A], "_vs_", levels(cluster_metadata$group_id)[B], "_all_"), quote = FALSE, = FALSE) #
Subset the significant results sig_res <- dplyr::filter(res_tbl, padj < padj_cutoff) %>%
dplyr::arrange(padj) (sig_res, paste0("DESeq2/pairwise/", clusters[x], "_", levels(cluster_metadata$group_id)[A], "_vs_", levels(cluster_metadata$group_id)[B], "_sig_"), quote = FALSE, = FALSE) ## ggplot of top genes normalized_counts <- counts(dds, normalized
= TRUE) ## Order results by padj values top20_sig_genes <- sig_res %>%
dplyr::arrange(padj) %>% dplyr::pull(gene) %>% head(n=20)
top20_sig_norm <- (normalized_counts) %>% rownames_to_column(var = "gene") %>% dplyr::filter(gene %in% top20_sig_genes) gathered_top20_sig <- top20_sig_norm %>% gather(colnames(top20_sig_norm)[2:length(colnames(top20_sig_norm))], key = "samplename", value = "normalized_counts")
gathered_top20_sig <- inner_join(ei[, c("sample_id", "group_id" )], gathered_top20_sig, by = c("sample_id" = "samplename")) ## plot using ggplot2 ggplot(gathered_top20_sig) + geom_point(aes(x = gene, y = normalized_counts,
color = group_id), position=position_jitter(w=0.1,h=0)) +
scale_y_log10() + xlab("Genes") + ylab("log10 Normalized Counts") +
ggtitle("Top 20 Significant DE Genes") + theme_bw() + theme(.x = element_text(angle = 45, hjust = 1)) + theme( = element_text(hjust = 0.5)) ggsave(paste0("DESeq2/pairwise/", clusters[x], "_", levels(cluster_metadata$group_id)[A], "_vs_", levels(cluster_metadata$group_id)[B], "_top20_DE_")) }# Run the script on all clusters comparing stim condition
relative to control conditionmap(1:length(clusters), get_dds_resultsAvsB, A = 2, B = 1)Script
to run DESeq2 on all cell type clusters - Likelihood Ratio TestThe following script will run th
版权声明：本文标题：DESeq2差异表达分析(二) 内容由网友自发贡献，该文观点仅代表作者本人， 转载请联系作者并注明出处：http://www.roclinux.cn/p/1703324530a446963.html， 本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容，一经查实，本站将立刻删除。
发表评论