转录组入门（mac 版本）

软件安装

  1. 安装bioconda:

  • 去官网下载和自己电脑系统一样的版本

https://conda.io/miniconda.html

  • 下载完后，双击解压，然后cd 到文件目录，开始安装。

  1. # 安装

  2. bash Miniconda3-latest-MacOSX-x86_64.sh

之后会要求看一个license，和安装环境路径，输入yes 就好。

  • 激活

  1. source .bash_profile

  2. # 查看是否安装成功

  3. conda

出现上图代表安装成功

  • bioconda 是conda 上一个分发生物信息的频道，所以需要创建一些channels.

  1. conda config --add channels conda-forge

  2. conda config --add channels defaults

  3. conda config --add channels r

  4. conda config --add channels bioconda

  • 安装完channels后 需要更新miniconda

  1. conda update conda

  1. 安装配置软件：

-sratoolkit

教程上用的Linux方法，mac会报错，建议使用以下方法。

方法一（别人可以，我的不行）：

  1. conda install -c bioconda sra-tools=2.8.1

方法二（亲测可行）：

  1. ruby -e “$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)”

  2. brew install wget

  3. brew tap homebrew/science

  4. brew install sratoolkit

  • 余下的安装很简单没报错

  1. conda install fastqc

  2. conda install hisat2

  3. conda install samtools

  4. conda install -c bioconda htseq

  5. conda install r

  6. conda install rstudio

参考文献

  1. http://www.biotrainee.com/thread-1796-1-1.html

  2. http://www.biotrainee.com/thread-1800-1-1.html

读文献下数据

  • 在NCBI查文章：

AKAP95 regulates splicing through scaffolding RNAs and RNA processing factors. Nat Commun 2016 Nov 8;7:13347. PMID: 27824034

  • 找数据地址

  • 去GSE网站，搜GSE81916

  • 下载代码

  1. # 在终端运行

  2. #用axel下载会快一些

  3. brew install axel

  4. for i in {56..62}

  5. do

  6. # 也可以用wget 下载

  7. axel ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP075/SRP075747/SRR35899$i/SRR35899$i.sra

  8. done

参考文献：

1. http://www.biotrainee.com/thread-1908-1-1.html

了解fastqc

需要用安装好的sratoolkit把sra文件转换为fastq格式的测序文件，并且用fastqc软件测试测序文件的质量！ 作业，理解测序reads，GC含量，质量值，接头，index，fastqc的全部报告，搜索中文教程，并发在论坛上面。 来源于生信技能树：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost

实验步骤

1. 将 sra 数据转化成 fastq 格式

先建立一个SRR_fastqc.sh 文件，写入代码

  1. #!/usr/bin/env bash

  2. for i in {56..62}

  3. do

  4. fastq-dump --gzip --split-3 -O /Users/chengkai/Desktop/zhuanlu_files -A SRR35899${i}.sra

  5. done

  • --gzip 压缩格式为gzip

  • --split-3 如果是双端测序输出两个文件，如果不是只输出一个文件

  • -0 输出文件路径

  • “/Users/chengkai/Desktop/zhuanlu_files” 这里改成你自己的文件路径

  • -A 输入文件路径

2. 在终端运行

  1. # 这一步很慢，我跑了2小时，泡杯咖啡，欣赏一部电影吧

  2. $ bash SRR_fastqc.sh 

3. fastqc 检测测序文件质量

创建一个文件夹

  1. mkdir fastqc/

创建一个fastqc.sh脚本，写入如下代码

  1. #!/usr/bin/env bash

  2. fastqc -o ./fastqc/ -t 8 SRR3589956_1.fastq.gz SRR3589956_2.fastq.gz

  3. fastqc -o ./fastqc/ -t 8 SRR3589957_1.fastq.gz SRR3589957_2.fastq.gz

  4. fastqc -o ./fastqc/ -t 8 SRR3589958_1.fastq.gz SRR3589958_2.fastq.gz

  5. fastqc -o ./fastqc/ -t 8 SRR3589959_1.fastq.gz SRR3589959_2.fastq.gz

  6. fastqc -o ./fastqc/ -t 8 SRR3589960_1.fastq.gz SRR3589960_2.fastq.gz

  7. fastqc -o ./fastqc/ -t 8 SRR3589961_1.fastq.gz SRR3589961_2.fastq.gz

  8. fastqc -o ./fastqc/ -t 8 SRR3589962_1.fastq.gz SRR3589962_2.fastq.gz

  • -o: 输出路径

  • -t: 线程数，和电脑配置有关，每个线程需要250MB 的内存

  1. bash fastqc.sh

4. 质量解读

html 格式用浏览器打开

基本信息

  • Enconding: 测序平台版本号

  • Total Sequence: reads 的数量

  • Sequence length: 总的序列数

  • %GC GC比，这个指标有物种意义，用于区别物种，一般人类42%

每个read各位置碱基的测序质量

  • 横轴碱基的位置（1-51），纵轴是质量分数，20表示1%的错误率，30表示0.1%

  • 红色线代表中位数，蓝色代表平均数，黄色是25%-75%区间，触须是10%-90%区间

  • Warning 报警 如果任何碱基质量低于10,或者是任何中位数低于25

  • Failure 报错 如果任何碱基质量低于5,或者是任何中位数低于20

偏离度

- 横轴碱基的位置（1-51） - 纵轴是tail的Index编号 - 检查reads中每一个碱基位置在不同的测序小孔之间的偏离度，蓝色代表偏离度小，质量好，越红代表偏离度越大，质量越差。 - 这个图主要是为了防止，在测序过程中，某些tail受到不可控因素的影响而出现测序质量偏低

reads质量的分布

  • 横轴表示Q值，0-40

  • 纵轴是每个值对应的reads数目

  • 当峰值小于27时，警告；当峰值小于20时，fail。我的报告峰值在38

GC 含量统计

- 横轴碱基的位置（1-51） - 纵轴是碱基含量百分比 - 图中四条线代表A T C G在每个位置平均含量 - 当部分位置碱基的比例出现bias时，即四条线在某些位置纷乱交织，往往提示我们有overrepresented sequence的污染。 - 本结果前10个位置，每种碱基频率有明显的差别，说明有污染。 - 当任一位置的A/T比例与G/C比例相差超过10%，报"WARN"；当任一位置的A/T比例与G/C比例相差超过20%，报"FAIL"

序列平均GC含量分布图

- 横轴是百分比； 纵轴是每条序列GC含量对应的数量

- 蓝色的线是程序根据经验分布给出的理论值，红色是真实值，两个应该比较接近才比较好 - 当红色的线出现双峰，基本肯定是混入了其他物种的DNA序列 - 偏离理论分布的reads超过15%时，报"WARN"；偏离理论分布的reads超过30%时，报"FAIL"

各位置N的reads比率

- 当测序仪器不能辨别某条reads的某个位置到底是什么碱基时，就会产生“N”，统计N的比率 - 正常情况下，N值非常小 - 当任意位置的N的比例超过5%，报"WARN"；当任意位置的N的比例超过20%，报"FAIL"

reads 长度分布

  • 每次测序仪测出来的长度在理论上应该是完全相等的

  • 当reads长度不一致时报"WARN"；当有长度为0的read时报“FAIL”

  • 当测序的长度不同时，如果很严重，则表明测序仪在此次测序过程中产生的数据不可信

统计不同拷贝数的reads的频率

- 横坐标是duplication的次数，纵坐标是duplicated reads的数目,以unique reads的总数作为100% - 测序深度越高，越容易产生一定程度的duplication，这是正常的现象，但如果duplication的程度很高，就提示我们可能有bias的存在

- 当非unique的reads占总数的比例大于20%时，报"WARN"；当非unique的reads占总数的比例大于50%时，报"FAIL"

接头含量

- 此图衡量的是序列中两端adapter的情况 - 如果在当时fastqc分析的时候-a选项没有内容，则默认使用图例中的四种通用adapter序列进行统计 - 本例中adapter都已经去除，如果有adapter序列没有去除干净的情况，在后续分析的时候需要先使用cutadapt软件进行去接头

重复短序列

  • 这个图统计的是，在序列中某些特征的短序列重复出现的次数

  • 我们可以看到1-8bp的时候图例中的几种短序列都出现了非常多的次数，一般来说，出现这种情况，要么是adapter没有去除干净，而又没有使用-a参数；要么就是序列本身可能重复度比较高，如建库PCR的时候出现了bias

  • 对于这种情况，我的办法是可以cut掉前面的一些长度，可以试着cut 1bp

参考文献

1. http://fbb84b26.wiz03.com/share/s/3XK4IC0cm4CL22pU-r1HPcQQ2irG2836uQYm2iZAyh1Zwf3_ （青山屋主） 2. www.biotrainee.com/thread-2034-1-1.html (laofuzi) 3. http://www.jianshu.com/p/14fd4de54402 （lxmic） 4. https://zhuanlan.zhihu.com/p/20731723 (孟浩巍)

了解参考基因和注释

在UCSC下载hg19参考基因组，从gencode数据库下载基因注释文件，并且用IGV去查看你感兴趣的基因的结构，比如TP53,KRAS,EGFR等等。 作业，截图几个基因的IGV可视化结构！还可以下载ENSEMBL，NCBI的gtf，也导入IGV看看，截图基因结构。了解IGV常识。

转录组文章分析思路

  • 好的文章不是靠写出来的，而是设计出来的。没有开始实验之前，就要先用一系列的假设去推导可能出现的结论，进行整体的构思。

  • 转录组测序的文章的核心内容是实验组与对照组的差异基因分析，无论是常见的发育调控、环境适应、免疫互作都离不开这部分核心的分析。

  1. 研究结果第一部分：转录组测序的实验设计中一般都会有一部分生理生化的实验结果，可能是电镜图片，可能是生理指标的变化，可能是有效物质的含量。

  2. 研究结果第二部分：根据生理生化指标选择了几个样本进行转录组测序？测了4G,6G数据量? 测序质量是否良好？检测到3000个基因表达？样本特异表达的基因有多少？

  3. 研究结果第三部分：差异分析组合间有多少差异基因？不同差异组合共有和特有的基因数量有多少？不同样本表达模式的差异怎样？

  4. 研究结果第四部分：差异分析组合间差异基因主要是那些基因家族？转录因子？KEGG通路？GO分类？这些研究结果需要去解释这生理生化指标的变化的原因

下载参考基因组

  • 登陆网站 http://genome.ucsc.edu/index.html 下载hg19参考基因组

  • hg下载前是900M 左右，解压后是3G大小

  1. # 下载文件

  2. axel http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz

  3. #  解压文件

  4. tar -zxvf chromFa.tar.gz

  5. #解压后可以发现，参考序列是按照染色体号分开列出的，我们还需要把所有的序列写入到一个文件中。

  6. cat *.fa > hg19.fa

  7. #最后删除其他无用的文件

  8. rm chr*.fa

下载注释文件

  • 官网：http://www.gencodegenes.org

  1. # GTF格式主要是用来描述基因的注释

  2. axel ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/GRCh37_mapping/gencode.v27lift37.annotation.gtf.gz

  3. # GFF文件是一种用来描述基因组特征的文件，现在我们所使用的大部分都是第三版）（GFF3）

  4. axel ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/GRCh37_mapping/gencode.v27lift37.annotation.gff3.gz

  5. # 解压并删除原来的文件

  6. gzip -d gencode.v27lift37.annotation.gtf.gz

  7. gzip -d gencode.v27lift37.annotation.gff3.gz

IGV软件的下载和安装

  • IGV软件全称：Intergrative Genomics Viewer 是一个高效的查看基因数据的可视化软件。

  • 官网：http://software.broadinstitute.org/software/igv/home

  • 安装igv 之前需要先安装java 8 以上版本

  • 官网：https://java.com/en/download/mac_download.jsp

IGV软件的使用

窗口

主窗口布局：

  1. tool bar（工具栏），menu bar（菜单栏），pop-up menus（弹出式菜单）

  2. 染色体上的红色盒子表示显示这部分染色体，显示完整染色体是红框会消失

  3. 尺度显示了染色体的可见部分，刻度线显示了染色体的位置，跨度列表显示了当前显示的碱基的数量

  4. IGV在水平行显示的数据称为tracks。通常,每个tracks代表一个样本或实验。这个例子展示了甲基化、基因表达、拷贝数，LOH和突变数据

  5. IGV也显示某些特性,比如在tracks中的基因。默认情况下,IGV在一个面板显示数据，在另一个面板显示数据特性。拖放一个track名称，将一个track从一个面板移动到另一个地方

  6. Track名称列在最左边面板。名字的易读性取决于 tracks的高度，例如，track越小，它的名字的可读性越小

  7. 属性名称被列在顶部的属性面板。彩色块代表属性值,每个独特的值被都有一个独特的颜色。鼠标放在一个颜色块的附近来查看其属性值

  8. 8.

导入参考基因组及注释信息，查看感兴趣基因的结构

  • 导入前面的 hg19.fa

  • 接着导入gtf文件，需要先sort，才能导入

-现在可以导入sort 后的gtf 文件了。

参考文献

  1. http://www.jianshu.com/p/48b5a0972301 (GTF/GFF文件的差异及其相互转换)

  2. http://www.jianshu.com/p/3e545b9a3c68 （hoptop）

  3. http://fbb84b26.wiz03.com/share/s/3XK4IC0cm4CL22pU-r1HPcQQ1lZQRc2nKQhn2SthRW24I8CZ （greenhillman）

序列比对

  • 比对软件很多，首先大家去收集一下，因为我们是带大家入门，请统一用hisat2，并且搞懂它的用法。 直接去hisat2的主页下载index文件即可，然后把fastq格式的reads比对上去得到sam文件。 接着用samtools把它转为bam文件，并且排序(注意N和P两种排序区别)索引好，载入IGV，再截图几个基因看看！ 顺便对bam文件进行简单QC，参考直播我的基因组系列。

下载 index

  1. axel ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz

  2. tar -zxvf hg19.tar.gz  #解压

  • ps.下了我19个小时。。。

正式比对

  • 命令行目录：Desktop/zhuanlu_files

  • hg19目录：Desktop/zhuanlu_files／hg19

  • RNA-Seq Data: Desktop/zhuanlu_files/RNA-Seq

  • 注释： -t 记录时间 -x hg19(index)文件路径 -1 -2 测序的两个fastq文件 -S 比对结果输出路径

  1. mkdir -p RNA-Seq/aligned

  2. for i in `seq 56 58`

  3. do

  4.    hisat2 -t -x hg19/genome -1 RNA-Seq/SRR35899${i}.sra_1.fastq.gz -2 RNA-Seq/SRR35899${i}.sra_2.fastq.gz -S RNA-Seq/aligned/SRR35899${i}.sam

  5. done

  • i7处理器，内存16G 运行结果如下，如果你在服务器上跑，内存大于64G，上面代码.sam后加一个'&’ 那就是6个线程同时跑了，大大提升速度，

  • 比对解读

    • 总共比对了28856780条，7.4%一次都没有比对到，77.18%比对到一次，15.41%大于一次

    • 没有匹配上的7.4%，不按照顺序再匹配，有4.3%匹配了

    • 之后把剩余部分用单端数据进行比对，64.28%没有比对上，26.16比对了一次，9.56%大于一次

    • 总共比对到95.45%

SAMtools

SAMtools 详解： http://www.jianshu.com/p/15f3499a6469 SAM 数据格式是目前高通量测序中存放比对数据的标准格式，但是SAM 文件都很大，非常占空间，所以需要转到bam文件，而用的就是SAMtools软件。

  • 注释 -S 输入sam文件 -b 指定输出的文件为bam

  • 排序默认是根据染色体的位置

  1. for i in `seq 56 58`

  2. do

  3.    samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam

  4.    samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam

  5.    samtools index SRR35899${i}_sorted.bam

  6. done

IGV 可视化

  • 文件导入“详见转录组入门四”

reads 计数

  1. 实现这个功能的软件也很多，还是烦请大家先自己搜索几个教程，入门请统一用htseq-count，对每个样本都会输出一个表达量文件。

  2. 需要用脚本合并所有的样本为表达矩阵。参考：生信编程直播第四题：多个同样的行列式文件合并起来。

  3. 对这个表达矩阵可以自己简单在excel或者R里面摸索，求平均值，方差。看看一些生物学意义特殊的基因表现如何，比如GAPDH,β-ACTIN等等。

分析前解读

原文中人类293cell的数据只有3个（SRR3589956-58）, 缺少一个对照重复。所以后续我只分析小鼠的四个样品（SRR3589959-62）

理论基础

  • 因表达定量有以下三个水平：

  • 基因水平

  • 转录水平

  • 外显子水平

  • 标准化

  • reads计数后得到的基因定量结果（count matrix），在进行不同维度的比较时需要进行不同的处理

  • 比较同一个样本(within-sample)不同基因之间的表达情况---主要考虑转录本长度的影响

  • 比较不同样本（across-sample）同一个基因的表达情况---主要考虑测序深度的影响

  • 标准化的算法：RPKM（SE）, FPKM（PE），TPM, TMM，RSEM等等。

  • 标准化之后才能进行差异表达分析

1. 下载小鼠index 文件和基因组注释文件

  1. axel ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz

  2. axel ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_mouse/release_M10/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf.gz

2. 序列与index比对

  1. for i in `seq 59 62`

  2. do

  3.    hisat2 -t -x mouse/mm10/genome -1 RNA-Seq/SRR35899${i}.sra_1.fastq.gz -2 RNA-Seq/SRR35899${i}.sra_2.fastq.gz -S RNA-Seq/aligned/SRR35899${i}.sam

  4. done

3. 用reads 排序bam 文件

  • 先将sam 文件转成bam文件，bam 文件用 read name (-n)排序，再转回到sam文件。

  1. for i in `seq 59 62`

  2. do

  3.    samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam

  4.    samtools sort -n SRR35899${i}.bam -o SRR35899${i}_nsorted.bam

  5.    samtools view -h SRR35899${i}_nsorted.bam > SRR35899${i}_count.sam

  6. done

4. HTSeq 进行基因水平的定量

HTSeq 详解

  • htseq-count 是一款用于reads计数的软件，他能对位于基因组上的一些单位的reads数进行统计，这里所说的单位主要是指染色体上的一组位置区间（我们常见的就是gene exon）

基本用法

  1. htseq-count [options] <alignment_file> <gff_file>

  1. :   比对到基因组后得到的SAM文件 （SAMtools包含一些perl 脚本可以将大多数的比对文件转换成SAM格式 ），注意基因组mapping时一定要用支持剪接的比对软件（splicing-aware aligner）进行比对软件如TopHat.  HTSeq-count 需要用到SAM格式中的 CIGAR 区域的信息。

  2. 想要通过标准输入来传入 基因组mapping得到SAM文件，用 – 替换  即可

  3. 如果你是双端测序，必须要对SAM进行排序（单端可不必排序，但这里我也推荐对单端测序结果排序已减少内存消耗并提高软件效率），对read name或 染色体位置 进行排序皆可（这里我推荐按read name排序，因为通过位置排序我遇到过错误）。具体需要通过 -r 参数指定，所以排序请详见参数 -r

  4. :  包含单位信息的gff/GTF文件（gff文件格式），大多数情况下就是指注释文件;   由于GTF文件其实就是gff文件格式的变形，在这里同样可以传入GTF格式文件。

模型

  • 如何判断一个reads属于某个基因， htseq-count 提供了union, intersectionstrict, intersectionnonempty 3 种模型，如图（大多数情况下作者推荐用union模型）：

参数

  • 指定输入文件的格式， 可以是 sam (for text SAM files) 或 bam (for binary BAM files)格式，默认是sam.

  1. -f <format>, –format=<format> ：

  • 对于双端测序数据，必须要对SAM文件进行排序，对read name或 位置 进行排序皆可，通过 -r 可以指定您的数据是以什么方式排序的： 可以是 name 或 pos ， 默认是name. 如果您的数据没有排序，可以通过samtools sort(推荐)或msort 排序， 命令：

  1. -r <order>, –order=<order>

  • 如果你指定name（按read name 排序）, htseq-count 期望输入的文件中的read pair是紧邻的2行；而指定pos (按位置排序) ，则不需要这样，进一步说，如果指定pos , 可以不排序，但这样会耗费更多的内存和效率。

  1. samtools sort -n accepted_hits_unique.bam -o accepted_hits_unique_name_sorted.bam   # -n 按read name 排序 ，如果不指定则按染色体位置排序

  • 是否这个数据是来自链特异性建库（默认 yes)

  1. -s <yes/no/reverse>, –stranded=<yes/no/reverse>

For stranded=no, a read is considered overlapping with a feature regardless of whether it is mapped to the same or the opposite strand as the feature. For stranded=yesand single-end reads, the read has to be mapped to the same strand as the feature. For paired-end reads, the first read has to be on the same strand and the second read on the opposite strand. For stranded=reverse, these rules are reversed.

  • 指定一个最低 read mapping质量值，低于值会被过滤掉（默认是10，0.5.4版本以前默认值是

  1. -a <minaqual>, –a=<minaqual>

  • 指定最小计数单位类型（gff文件的第3列中的类型如： exon）, 指定后其他单位类型将被忽略（默认情况下，对于rna-seq分析采用 Ensembl GTF 文件类型时,默认值是：exon）。

  1. -t <feature type>, --type<feature type>

  • GFF文件的一类属性，最终的技术单位，默认情况下，对于rna-seq分析采用 Ensembl GTF 文件类型时,默认值是：gene_id）。

  1. -i <id attribute>, --idattr=<id attribute>

  • 判断一个reads属于某个基因的模型，用来判断统计reads的时候对一些比较特殊的reads定义是否计入。 包括：默认的union和intersection-strict、 intersection-nonempty （默认：union）。

  1. -m <mode>, --mode=<mode>

  • 输出所有alignment的reads到一个叫 的sam文件中，通过一个可选的sam标签 ' XF ’来标注每一行对应的单位和计数，可以不设置。

  1. -o <samout>, --samout=<samout>

  • 屏蔽程序报告和警告

  1. -q, --quiet

  • 打印帮助

  1. -h, --help

5. reads 计数

  • 数据准备已经完成，接下来要使用htseq-count 对数据进行reads计数。

  1. # 安装HTSeq

  2. conda install HTSeq

  3. # 语法

  4. Usage: htseq-count [options] <sam_file> <gff_file>

  1. mkdir -p RNA-Seq/matrix/

  2. htseq-count RNA-Seq/aligned/SRR3589959_count.sam mouse/gencode.vM10.annotation.gtf > RNA-Seq/matrix/SRR3589959.count

  3. htseq-count RNA-Seq/aligned/SRR3589960_count.sam mouse/gencode.vM10.annotation.gtf > RNA-Seq/matrix/SRR3589960.count

  4. htseq-count RNA-Seq/aligned/SRR3589961_count.sam mouse/gencode.vM10.annotation.gtf > RNA-Seq/matrix/SRR3589961.count

  5. htseq-count RNA-Seq/aligned/SRR3589962_count.sam mouse/gencode.vM10.annotation.gtf > RNA-Seq/matrix/SRR3589962.count

  • 我是分四个终端一起跑，相当于四个进程同时跑

  • 最后生成四个文件

  • count文件的格式：基因名一列（ENSMUSG00000000001.4）+reads计数一列(1648)

6. 合并表达矩阵

  • 使用R将这四个文件进行合并（59 和 61 是对照， 60和62是处理），得到最后的融合表达矩阵，R语言代码如下

  1. options(stringsAsFactors = FALSE)

  2. # 首先将四个文件分别赋值：control1，control2，rep1，rep2

  3. control1 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589959.count", sep="\t", col.names = c("gene_id","control1"))

  4. control2 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589961.count", sep="\t", col.names = c("gene_id","control2"))

  5. rep1 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589960.count", sep="\t", col.names = c("gene_id","akap951"))

  6. rep2 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589962.count", sep="\t",col.names = c("gene_id","akap952"))

  7. # 将四个矩阵按照gene_id进行合并，并赋值给raw_count

  8. raw_count <- merge(merge(control1, control2, by="gene_id"), merge(rep1,rep2, by="gene_id"))

  9. # 需要将合并的raw_count进行过滤处理，里面有5行需要删除的行，在我们的小鼠的表达矩阵里面，是1,2,48823,48824,48825这5行。并重新赋值给raw_count_filter

  10. raw_count_filt <- raw_count[-48823:-48825, ]

  11. raw_count_filter <- raw_count_filt[-1:-2, ]

  12. # 因为我们无法在EBI数据库上直接搜索找到ENSMUSG00000024045.5这样的基因，只能是ENSMUSG00000024045的整数，没有小数点，所以需要进一步替换为整数的形式。

  13. # 第一步将匹配到的.以及后面的数字连续匹配并替换为空，并赋值给ENSEMBL

  14. ENSEMBL <- gsub("\\.\\d*", "", raw_count_filter$gene_id)

  15. # 将ENSEMBL重新添加到raw_count_filter矩阵

  16. row.names(raw_count_filter) <- ENSEMBL

  17. # 看一些基因的表达情况，在UniProt数据库找到AKAP95的id，并从矩阵中找到访问，并赋值给AKAP95变量

  18. AKAP95 <- raw_count_filter[rownames(raw_count_filter)=="ENSMUSG00000024045",]

  19. # 查看AKAP95

  20. AKAP95

7. 简单分析

  • 合并后的数据框

  • AKAP95 处理组明显比对照组高，说明这个基因在小鼠的表达量变高

  • 总的概况

8. 参考文献

  1. http://www.jianshu.com/p/e9742bbf83b9 （hoptop 大神）

  2. http://www.jianshu.com/p/24cf44b610a7 （lxmic 大神）

  3. http://fbb84b26.wiz03.com/share/s/3XK4IC0cm4CL22pU-r1HPcQQ3ouv263iAk012b-9tU2UXRgt（青山屋主大神）

差异基因分析

  • 差异基因分析这个步骤推荐在R里面做，载入表达矩阵，然后设置好分组信息，统一用DEseq2进行差异分析，当然也可以走走edgeR或者limma的voom流程。基本任务是得到差异分析结果，进阶任务是比较多个差异分析结果的异同点。生信技能树

1. 读取表达矩阵

  1. # 首先将四个文件分别赋值：control1，control2，rep1，rep2

  2. control1 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589959.count", sep="\t", col.names = c("gene_id","control1"))

  3. control2 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589961.count", sep="\t", col.names = c("gene_id","control2"))

  4. rep1 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589960.count", sep="\t", col.names = c("gene_id","akap951"))

  5. rep2 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589962.count", sep="\t",col.names = c("gene_id","akap952"))

  6. # 将四个矩阵按照gene_id进行合并，并赋值给raw_count

  7. raw_count <- merge(merge(control1, control2, by="gene_id"), merge(rep1,rep2, by="gene_id"))

  8. # 需要将合并的raw_count进行过滤处理，里面有5行需要删除的行，在我们的小鼠的表达矩阵里面，是1,2,48823,48824,48825这5行。并重新赋值给raw_count_filter

  9. raw_count_filt <- raw_count[-48823:-48825, ]

  10. raw_count_filter <- raw_count_filt[-1:-2, ]

  11. # 因为我们无法在EBI数据库上直接搜索找到ENSMUSG00000024045.5这样的基因，只能是ENSMUSG00000024045的整数，没有小数点，所以需要进一步替换为整数的形式。

  12. # 第一步将匹配到的.以及后面的数字连续匹配并替换为空，并赋值给ENSEMBL

  13. ENSEMBL <- gsub("\\.\\d*", "", raw_count_filter$gene_id)

  14. row.names(raw_count_filter) <- ENSEMBL

  15. raw_count_filter <- raw_count_filter[ ,-1]

2. DESeq2

  1. # 下载DESeq2

  2. source("https://bioconductor.org/biocLite.R")

  3. biocLite("GenomeInfoDb")

  • DESeq2 对于输入数据的要求

  1. DEseq2要求输入数据是由整数组成的矩阵。

  2. DESeq2要求矩阵是没有标准化的。

  • DESeq2进行差异表达分析

  • DESeq2包分析差异表达基因简单来说只有三步：构建dds矩阵，标准化，以及进行差异分析。

  1. dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition) #~在R里面用于构建公式对象，~左边为因变量，右边为自变量。

  2. dds <- DESeq(dds) #标准化

  3. res <- results(dds, contrast=c("condition","treated","control")) #差异分析结果

3. 构建dds矩阵

  • 表达矩阵：即上述代码中的countData，就是我们前面通过read count计算后并融合生成的矩阵，行为各个基因，列为各个样品，中间为计算reads或者fragment得到的整数。

  • 样品信息矩阵：即上述代码中的colData，它的类型是一个dataframe（数据框），第一列是样品名称，第二列是样品的处理情况（对照还是处理等），即condition，condition的类型是一个factor。

  • 差异比较矩阵： 即上述代码中的design。 差异比较矩阵就是告诉差异分析函数是要从要分析哪些变量间的差异，简单说就是说明哪些是对照哪些是处理。

  1. # 这一步很关键，要明白condition这里是因子，不是样本名称；小鼠数据有对照组和处理组，各两个重复

  2. > condition <- factor(c(rep("control",2),rep("akap95",2)), levels = c("control","akap95"))

  3. # 获取count数据

  4. > countData <- raw_count_filter[,1:4]

  5. > colData <- data.frame(row.names=colnames(raw_count_filter), condition)

  6. > dds <- DESeqDataSetFromMatrix(countData, colData, design= ~ condition)

  7. # 查看一下dds的内容

  8. > head(dds)

4. DESeq标准化dds

  1. # normalize 数据

  2. > dds2 <- DESeq(dds)

  3. # 查看结果的名称，本次实验中是 "Intercept"，"condition_akap95_vs_control"

  4. > resultsNames(dds2)

  5. # 将结果用results()函数来获取，赋值给res变量

  6. res <- results(dds2)

  7. # summary一下，看一下结果的概要信息

  8. summary(res)

5. 提取差异分析结果

  1. # 获取padj（p值经过多重校验校正后的值）小于0.05，表达倍数取以2为对数后大于1或者小于-1的差异表达基因。

  2. > table(res$padj<0.05)

  3. > res <- res[order(res$padj),]

  4. > diff_gene_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))

  5. > diff_gene_deseq2 <- row.names(diff_gene_deseq2)

  6. > resdata <-  merge(as.data.frame(res),as.data.frame(counts(dds2,normalize=TRUE)),by="row.names",sort=FALSE)

  7. # 得到csv格式的差异表达分析结果

  8. > write.csv(resdata,file= "control_vs_akap95.cvs",row.names = F)

6. 参考文献

  1. https://zhuanlan.zhihu.com/p/30350531（青山屋主）

  2. http://www.jianshu.com/p/3bfb21d24b74 （lxmic）

  3. http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html （DESeq2 官网）

赞 (0)

相关推荐