Category Archives: 未分类

六 25

这是一个违背祖宗的决定

Posted on by

最近看到朋友圈转发的一大批“神医”做出违背祖宗的决定！

腾讯视频链接：https://v.qq.com/x/page/x3230xgj0x6.html

让我动容，敬佩之外我也想效仿一二，把我“祖传的”生物信息学技能公之于众，敞开门让大家学！我已经把这些技能录制成为了视频，并且在B站免费发布，已经组建了微信交流群的有下面这些：

当然了，如果你直接上手这些NGS组学数据分析实战有困难，说明你的Linux基础或者R语言不过关，也可以看c超级基础的内容：

image-20210303182603072

再怎么强调生物信息学数据分析学习过程的计算机基础知识的打磨都不为过，我把它粗略的分成基于R语言的统计可视化，以及基于Linux的NGS数据处理：

把R的知识点路线图搞定，如下：

  • 了解常量和变量概念
  • 加减乘除等运算（计算器）
  • 多种数据类型（数值，字符，逻辑，因子）
  • 多种数据结构（向量，矩阵，数组，数据框，列表）
  • 文件读取和写出
  • 简单统计可视化
  • 无限量函数学习

Linux的6个阶段也跨越过去 ，一般来说，每个阶段都需要至少一天以上的学习：

  • 第1阶段：把linux系统玩得跟Windows或者MacOS那样的桌面操作系统一样顺畅，主要目的就是去可视化，熟悉黑白命令行界面，可以仅仅以键盘交互模式完成常规文件夹及文件管理工作。
  • 第2阶段：做到文本文件的表格化处理，类似于以键盘交互模式完成Excel表格的排序、计数、筛选、去冗余，查找，切割，替换，合并，补齐，熟练掌握awk,sed,grep这文本处理的三驾马车。
  • 第3阶段：元字符，通配符及shell中的各种扩展，从此linux操作不再神秘！
  • 第4阶段：高级目录管理：软硬链接，绝对路径和相对路径，环境变量。
  • 第5阶段：任务提交及批处理，脚本编写解放你的双手。
  • 第6阶段：软件安装及conda管理，让linux系统实用性放飞自我。

文末友情推荐

​
六 25

seq-ImmuCC实战推断小鼠的肿瘤免疫微环境

Posted on by

前面我在《生信菜鸟团》公众号介绍了 小鼠的肿瘤免疫微环境推断可以用seq-ImmuCC，本来是想布置作业让大家下载https://www.ncbi.nlm.nih.gov/bioproject/PRJNA489661/ 的数据，自己走一波seq-ImmuCC实战。
但是发现这个并不是《生信技能树》公众号，所以我没办法布置练习题额，所以现在在《生信技能树》公众号再来一次哈！前面的例子：人人都能学会的单细胞聚类分群注释 ，第一次分群就非常漂亮！可以看到这个数据集GSE129516里面的6个样品都是有不同的免疫细胞亚群的，而且既然已经是有了降维聚类分群结果，就可以算出各个细胞亚群的比例啦！
我们的作业是，把GSE129516里面的6个样品的单细胞表达量矩阵简单的累加成为一个假的bulk表达量矩阵， 然后拿这个表达量矩阵去进行seq-ImmuCC实战，推断小鼠的肿瘤免疫微环境，就是各个免疫细胞比例。然后跟单细胞的真实免疫细胞各个亚群比例进行对比，测试一下seq-ImmuCC这个网页工具的表现情况！

关于seq-ImmuCC开发者

苏州系统医学研究所苏吴爱平教授和秦晓峰教授合作，发布了一个基于转录组测序（RNA-Seq）数据对小鼠组织中10种主要免疫细胞组分进行预测的计算模型（seq-ImmuCC），将为测序数据提供免疫细胞层面的解读视角。相关研究以“seq-ImmuCC: Cell-CentricView of Tissue Transcriptome Measuring Cellular Compositions of ImmuneMicroenvironment From Mouse RNA-Seq Data”为题，于2018年6月发表于国际期刊Frontiersin Immunology。
吴爱平教授的课题组已经成功开发了两个计算模型，从组织样本的芯片表达谱。

​
六 25

RNA-Seq数据的内含子保留分析

Posted on by

前面我们的明码标价之普通转录组上游分析，受到了各大热心粉丝的吐槽，觉得太简单了我们居然还好意思收费。后面我就就加上了稍微有一点难度的《可变剪切》，不过仍然是阻挡不了粉丝无穷无尽的需求，后台有人发给我一个RNA-Seq数据的内含子保留分析需求。

我看了看粉丝发过来的文章，发表于 January 2021, 在CELL杂志的文章《Spliceosome-targeted therapies trigger an antiviral immune response in triple-negative breast cancer》，链接是：https://doi.org/10.1016/j.cell.2020.12.031

这个文章数据比较多：

SUM159 SD6 RNA-Seq #GSE163414
LM2 SD6 RNA-Seq #GSE163411
SUM159 Cytoplasmic RNA-Seq #GSE163232
SUM159 J2 dsRIPseq #GSE163188
Syngeneic model RNA-Seq #GSE163181

可以看到，主要是RNA-Seq数据啦，有两个是普通的细胞系处理前后的表达量差异情况探索，所以出图如下：

表达量差异情况探索标准图表

这个已经是超级简单了， 我们的明码标价之转录组常规测序服务（仅需799每个样品） 和 明码标价之普通转录组上游分析 就是对这样的 RNA-Seq拿到了表达量矩阵，然后下游分析也是平淡无奇，仅收费800，代码呢，我也多次分享了，基本上看我六年前的表达芯片的公共数据库挖掘系列推文即可；

这样的分析流程基本上绝大部分粉丝已经是无需委托我们啦，所以粉丝发给我的是 RNA-Seq数据的内含子保留分析需求，步骤如下：

  • Hisat2-aligned reads were filtered for proper-paired reads (-f 2 flag in SAMtools).
  • Intron annotations were parsed from UCSC RefSeq gene annotation files and were filtered to exclude features that overlap genomic loci on the same strand.
  • Reads mapping to introns were counted using Pysam.
  • For each intron feature, we defined the following two read classes:
    • (1) ‘‘intronic’’ reads mapping at least 6 bases contiguously within the intron
    • (2) ‘‘spanning’’ reads with ends mapping to the flanking exons.
  • The intron retention (IR) score was then computed as the ratio of the RPKM-normalized ‘‘intronic’’ read density over the RPKM-normalized ‘‘spanning’’ read den- sity.
  • In order to compare commonly expressed IR events across samples, introns with < 10 spanning RPKM in any sample were excluded from all analyses.

对我们有ngs组学数据分析经验的人来说，其实并不难，无非就是安装几个软件，使用几个包。但对于没有学过编程的纯生物学研究者来说基本上不可能完成，也没有这样的网页工具。

但是呢，这个流程又确实是过于个性化，哪怕对我们来说很简单，也其实是耗费时间和精力需要研发调试的。

首先你需要有RNA-seq的fastq文件

如果是TCGA数据库，步骤如下：

  • Intron retention analysis was performed on BRCA TCGA RNA sequencing datasets (Koboldt et al., 2012).
  • TCGA fastq reads were mapped using the STAR aligner (v2.3.1) (Dobin et al., 2013) onto the hg19/GRCh37 reference genome as previously described (Hsu et al., 2015).
  • Level of intron retention (IR level) within each sample was calculated as the number of introns with IR scores > 0.01, as defined previously.
  • ‘‘High’’ and ‘‘Low’’ IR were defined as having an IR level outside one standard deviation of the mean.
  • RSEM normalized gene expression data from TCGA was obtained from the Broad GDAC Firehose.

一般来说，大家是很难下载TCGA数据库原始fastq文件，这个权限审核比较严厉，不过咱们数据挖掘呢完全没有毕业一直盯着TCGA数据库啊，自己领域的普通RNA-seq肯定也是不少。如果是认真搞科研，你一定会自行调研和阅读文献，找到合适的数据集。

数据挖掘的核心就是通过分类把基因数量搞少

部分粉丝看到这里，可能无法理解RNA-Seq数据的内含子保留分析的意义是什么？

其实就是多了一个维度的指标，来把你的样本分类，分类后就可以找差异。同样的我们可以看这个示例文章，感觉每个样品的IR指标，把病人分成IR高低两个组别，然后走普通的ssGSEA分析，生存分析。

image-20210320093952534

这一套组合拳，大家是不是很眼熟啊？

如果你也想做自己的的RNA-Seq数据的内含子保留分析，赶快联系我们吧，同样的，我们的分析仍然是明码标价，单个RNA-Seq数据的内含子保留分析收费仅需800元，因为是纯粹的基于Linux平台的各种软件脚本，所以提供你全套数据和脚本但是无法保证你能运行成功，因为你不一定有自己的服务器。

​
六 25

去broad官网下载msigdb数据库文件很麻烦

Posted on by

我在：借鉴escape包的一些可视化GSVA或者ssGSEA结果矩阵的方法 和 对单细胞表达矩阵做gsea分析的两个教程里面提到过，MSigDB（Molecular Signatures Database）数据库中定义了已知的基因集合：http://software.broadinstitute.org/gsea/msigdb 需要注册才能下载。

但是这个GitHub包，ncborcherding/escape文档，在：http://www.bioconductor.org/packages/release/bioc/vignettes/escape/inst/doc/vignette.html 提供了一个封装好的MSigDB数据库信息，其实你仔细看它的文档，它的打包其实是依赖于msigdbr_7.2.1。

获取 MigDB中的全部基因集

MigDB中的全部基因集 都被这个GitHub包，ncborcherding/escape 打包起来了，MSigDB（Molecular Signatures Database）数据库中定义了已知的基因集合：http://software.broadinstitute.org/gsea/msigdb 包括H和C1-C7八个系列（Collection），每个系列分别是：

  • H: hallmark gene sets （癌症）特征基因集合，共50组，最常用；
  • C1: positional gene sets 位置基因集合，根据染色体位置，共326个，用的很少；
  • C2: curated gene sets：（专家）校验基因集合，基于通路、文献等：
  • C3: motif gene sets：模式基因集合，主要包括microRNA和转录因子靶基因两部分
  • C4: computational gene sets：计算基因集合，通过挖掘癌症相关芯片数据定义的基因集合；
  • C5: GO gene sets：Gene Ontology 基因本体论，包括BP（生物学过程biological process，细胞原件cellular component和分子功能molecular function三部分）
  • C6: oncogenic signatures：癌症特征基因集合，大部分来源于NCBI GEO 发表芯片数据
  • C7: immunologic signatures: 免疫相关基因集合。
GS <- getGeneSets(library = "H")
GS

MigDB中的全部基因集 被 构建成为： a list of GSEABase GeneSet objects ，获取 hallmark gene sets （癌症）特征基因集合。

源头是msigdbr 包

安装方法非常简单：

install.packages("msigdbr")

但是这个msigdbr并没有我想象中的那么大：

Installing package into ‘C:/Users/win10/Documents/R/win-library/4.0’
(as ‘lib’ is unspecified)
试开URL’https://cran.rstudio.com/bin/windows/contrib/4.0/msigdbr_7.2.1.zip'
Content type 'application/zip' length 6737651 bytes (6.4 MB)
downloaded 6.4 MB

package ‘msigdbr’ successfully unpacked and MD5 sums checked

同样的，学习R包，看看文档即可，在： https://cran.r-project.org/web/packages/msigdbr/vignettes/msigdbr-intro.html

Documentation for package ‘msigdbr’ version 7.2.1
DESCRIPTION file.
User guides, package vignettes and other documentation.
Help Pages
msigdbr Retrieve the gene sets data frame
msigdbr_collections List the collections available in the msigdbr package
msigdbr_show_species List the species available in the msigdbr package
msigdbr_species List the species available in the msigdbr package

非常简单的文档

这些代码使用就明白了，确实没啥好继续讲解的：

library(msigdbr)
# All gene sets in the database can be retrieved without specifying a collection/category. 
all_gene_sets = msigdbr(species = "Mus musculus")
head(all_gene_sets)
msigdbr_species()
all_gene_sets = msigdbr(species = "Homo sapiens")

无非就是封装和对象，前面的 escape 包提供了getGeneSets函数，我们的这个msigdbr提供了 msigdbr函数。

生信基石之R语言

B站的10个小时教学视频务必看完，参考 GitHub 仓库存放的相关学习路线指导资料：https://github.com/jmzeng1314/R_bilibili ，可以参考一些优秀笔记，比如https://mubu.com/doc/2KUiSCfVsg

文末友情推荐

​