RNA-seq实战（1）

本来想往下讲讲机器学习的，但是前段时间有一些事情耽搁了。然后我越学越发现机器学习深度之深，广度之广。包括在python中Tensorflow构建神经元等，我感觉目前以我的理解和我的能力，只能解决一些简单问题而不能解决一些系统问题，所以我还是将其先搁置。待我对整个machine learning 有一个更深入的认识后再和大家分享吧。

所以今天就讲讲最近练手分析的一波数据吧。我之后的所有内容将围绕这篇2017年发表在JCI insight上的杂志来展开

虽然，JCI insight是一个16新办的杂志，还没有影响因子，不过我看了一下引用率还是很不错的。所以这项研究还是不错的，作为我们练习的素材还是ok的。

这篇文章主要鉴定了SMOC2，一种胞外分泌钙调蛋白，在肾小管纤维化过程中的生物学作用

pro：大概是一个正相关关系，建立的SMOC2超表达小鼠，基本都出现了肾小管纤维化现象

con： 敲除SMOC2的细胞系基本无法诱导肾小管细胞纤维化，siRNA干扰SMOC2也无法诱导肾小管细胞纤维化

那么作者为什么会关注SMOC2这个基因呢。原因是他们前一年对纤维化肾小管细胞做RNA-seq测序时，发现SMOC2有很高的表达。所以将SMOC2拿出来研究。

我们直接切入RNA-seq的部分吧。为了研究SMOC2下游调控的基因。作者构建了超表达SMOC2的小鼠，根据其是否产生纤维化肾小管分为normal和fibrosis。

包括实验对照组在内，共有14个样本。

在RNA-seq实验设计中有几个需要考虑的重要因素。一个是生物学重复（biological replicate），和批次效应（batch effect）因为篇幅的原因我在这里不展开。请大家自行学习

在测序上机结束后，我们会得到很大的fastq文件。他大概长什么样子呢。

第一行代表了测序序列的名字和一些其他信息，第一行最后一个一是什么意思呢？是其测序不同链的记号。（RNA-seq测序一般会有postive strand forward，negative strand forward）第二行就是我们的测序核酸序列，可以看到序列中存在很多的NNNN，那个代表这里的测序结果并不好。如果第四行像我们的结果一样均为####则代表这里的测序结果真的很差。对fastq文件进行质量控制的很简单，在这里我也不展开了。

那么我们知道RNA-seq是将cDNA片段化（fragmented），再测序（二代测序基本原理不理解的话建议上Youtube搜一下相关视频，有很多很详细的的动画解答），然后再比对到基因组（当然也有不利用基因组直接比对的，但这不是我们的重点）

我们利用很不同的alignment package就可以将我们序列比对到基因组（因为基本pipeline基本相同在这里不想展开）。这时候我们就可以得到gene count table.

在文章中我们可以找到作者上传到GEO 数据库的序列号，在GEO数据中的均为可以直接进行下游分析的Gene count table

我们到GEO数据库找到我们的序列号。

可以看到他的简介和他的数据，我们可以看到一共上传了14个测序数据，即我们的十四个样本。

我们把它们手动下载下来，然后放到我们的工作环境中。

当然大家也可以用GEOquery packages来下载。不过目前对下载下来的矩阵如何提取结果我还遇到一些小问题尚未解决，所以我就先讲讲我自己开发的土方法吧。虽然真的很土，但是也还算好用（捂脸）。

我们主要用到两个包，一个stringr，一个dplyr

第一步我先将workspace中的所有文件变量名用dir（）提取。然后删除掉文件名中的冗余部分

其完整的文件名包涵样本信息和序列号等，我主要用stringr中的str_sub删除前面12个字符，后面15个字符就可以得到sample的名字

d

然后我们就可以得到这么一个gene count table，列上均为我们的实验条件（基因型+是否纤维化），而行均为ensemble上的基因accession number。然后另一个就是实验的metadata，什么是metadata，就是描述数据的数据（the data that describe the data）我们再来看看实验的条件

我们有七个野生型样本，七个SMOC2超表达样本。其中都是四个纤维化样本，四个正常样本

根据这个我来创建metadata。都是一些常规操作，如果有不太理解的欢迎看一下我之前写的关于R的基本文章

然后这就是我们创建的metadata，还有一步工作就是要把我们rawdata的顺序和metadata的顺序对上

我们可以用all（）来看看

如果他们的顺序不一样我们就要重新重排数据，有一个小细节就是第一步提取了geneID作为行名，然后删除了rawdata中的列geneID，这主要是为了让rawdata和metadata保持一样的列名，在我们分析完之后又要重新提取出来（后面讲）。

然后我们对数据的一些前期处理和准备就做好了。下一步就是要对数据进行质量控制（quality control），然后要用到层序聚类分析和PCA这两个非监管学习方法。我在之前写过层序聚类分析，如果有不太了解的同学，还要看一下前面的文章