NGS研究策略和应用技术-转录组上机文档

转录组分析上机文档

一、 过putty登录服务器

1. 通过SSH客户端putty登录服务器，以及安装xming

点击链接http://the.earth.li/~sgtatham/putty/latest/x86/putty.exe下载putty.exe，点击链接http://jaist.dl.sourceforge.net/project/xming/Xming/6.9.0.31/Xming-6-9-0-31-setup.exe下载xming，至任意位置，安装xming后运行程序。

双击打开putty，在Connection-SSH-X11栏，勾选“enable X11 forwarding”,在x11 display在location中填入“localhost:0”；在Session栏的主机地址输入服务器地址192.168.30.31，点击Load打开。如果出现安全警告框，点击YES。再输入用户名和密码登陆服务器。

设置linux的可视化界面：

登陆linux终端：

终端进入，欢迎界面：

二、 软件和数据

1) 分析用到的软件

基因组比对部分用到的软件是：

bowtie[官网：http://bowtie-bio.sourceforge.net/bowtie2/index.shtml] tophat[官网：http://tophat.cbcb.umd.edu/]； 表达量估算用到的软件是：

cufflinks[官网：http://cufflinks.cbcb.umd.edu/]；

差异分析用到的软件是cuffdiff,为cufflinks软件包的一部分。 分析软件均为开源软件，可以下载最新版本，进行安装。

2) 准备分析数据

在服务器的家目录下建立文件夹RNA-seq,然后建立软链接，把teacher/RNA-seq的数据，在本地建立快捷方式。此后的工作目录均为是在RNA-seq下进行。

命令行： mkdir RNA-seq ln -s /home/teacher/RNA-seq/data RNA-seq/ cd RNA-seq; ls -l data

打开data文件夹进行查看，

文件说明：

fq文件为测序数据，C460和P460表示测序样本条件1和测序样本条件2，R1.fq和R2.fq表示双末端pair-end数据；

GRCh37chr21.fa文件是human 染色体第21号染色体的fasta序列，用来做比对的参考基因组；

GRCh37chr21.fa.fai文件是基因组文件的索引文件。 genes.gtf文件是gtf格式的基因组结果注释文件；

三、 分析流程

1. 认识和查看FASTQ测序数据

利用基础的Linux命令，查看FASTQ格式的测序数据： head data/C460-R1.fq head 查看文件的前十行信息

FASTQ是基于文本的，保存生物序列（通常是核酸序列）和其测序质量信息的标准格式。其序列以及质量信息都是使用一个ASCII字符标示，最初由Sanger开发，目的是将FASTA序列与质量数据放到一起，目前已经成为高通量测序结果的事实标准。FASTQ文件中每个序列通常有四行：

1. 序列标识以及相关的描述信息，以‘@’开头； 2. 第二行是序列

3. 第三行以‘+’开头，后面是序列标示符、描述信息，或者什么也不加

4. 第四行，是质量信息，和第二行的序列相对应，每一个序列都有一个质量评分，根据评分体系的不同，每个字符的含义表示的数字也不相同。

2. 使用FastQC进行测序数据质量检测

cd$HOME/RNA-seq mkdirqc_out fastqcdata/*.fq -o qc_out 组装参数的说明： *.fq表示对当前目录下面的所有以fq为后缀名的文件进行质控，-o表示输出的目录。

结果：qc_out目录下会生成四个文件夹及四个压缩文件，四个文件夹的质控结果分别对应前面的四个fastq文件

3. 基因组比对

基因组比对是以参考基因组为模版，把测序的短片段reads比对到参考基因组的具体位置上。在转录组的分析中，我们使用tophat软件进行比对，由于tophat需要调用bowtie，在使用前需用bowtie建立参考基因组的索引文件。

i. 建立基因组索引文件 (约1min)

命令行： cd cd RNA-seq bowtie2-build data/GRCh37chr21.fa chr21 运行完后，在当前目录[RNA-seq目录]下进行查看，可看到生成以下结果文件：

参数说明：GRCh37chr21.fa 是基因组文件， chr21是基因组索引文件的前缀名称。 结果文件说明：

chr21.1.bt2; chr21.2.bt2; chr21.3.bt2; chr21.4.bt2; chr21.rev.1.bt2; chr21.rev.1.bt2; 是建立的索引文件，共6个。

ii. 用tophat进行基因组比对 (约12min/file)

命令行：

cd ~/RNA-seq/ 运行完后，在当前目录[RNA-seq目录]下生成文件夹C460_mapping 和P460_mapping。 tophat -p 1 -G data/genes.gtf -o C460_mapping chr21 data/C460-R1.fq data/C460-R2.fq tophat -p 1 -G data/genes.gtf -o P460_mapping chr21 data/P460-R1.fq data/P460-R2.fq 如进入C460_mapping，可以看到如下结果:

参数说明：-p 1 表示使用一个线程；-G genes.gtf表示使用已知基因结构注释文件，在比对中，bowtie会优先对这部分进行比对；-o 表示生成文件夹的名称; 其他参数全都使用默认参数。

结果文件说明：accepted_hits.bam表示比对上的reads生成的比对文件，unmapped.bam表

示没有比对上的reads生成的比对文件；这两个文件都可以用samtools工具进行查看[samtools view –h accepted_hits.bam|head -10]; deletions.bedinsertions.bedjunctions.bed是tophat对可变剪接位点的记录文件，分别表示没有检测到的可变剪接位点；识别出的剪接位点和所有已知未知的可变剪接位点,可以用less进行查看[ less junctions.bed|head -10 ]; align_summary.txt 是tophat2.0.09版本及以上对reads比对情况的统计；pre_reads.info和logs文件是对进程中的reads数目和进程运行情况的记录。

iii. 生成比对文件的索引文件(约1min)[可视化部分使用]

命令行：

cd ~/RNA-seq/ samtools index C460_mapping/accepted_hits.bam samtools index P460_mapping/accepted_hits.bam

运行完后，在当前目录[RNA-seq目录]的C460_mapping和P460_mapping文件夹下，分别生成文件.bam.bai。

如进入C460_mapping，可以看到如下结果,生成了文件accepted_hits.bam.bai:

4. 表达量估算 (约1min/file)

在表达量估算部分，我们使用cufflinks工具，先对短片段reads进行转录本重构，然后计算基因和转录本的表达量，以FPKM为单位。

命令行：

cd ~/RNA-seq/ cufflinks -p 1 -G data/genes.gtf -o C460_expression C460_mapping/accepted_hits.bam cufflinks -p 1 -G data/genes.gtf -o P460_expression P460_mapping/accepted_hits.bam

运行完后，在当前目录[RNA-seq目录]下，生成文件夹C460_expression和P460_expression。

如进入C460_expression文件夹，可以看到如下结果:

参数说明：-p 1 表示使用一个线程；-G genes.gtf表示使用已知基因结构注释文件，对这些有结构注视的基因和转录本进行表达量估算；-o 表示生成的文件夹名称。

结果文件说明：genes.fpkm_tracking是生成的基因表达量的文件；isoforms.fpkm_tracking是生成的转录本表达量的文件；transcripts.gtf, skipped.gtf分别表示检测了的基因结构文件和滤掉的基因结构文件。

5. 差异表达分析 (约2min)

在表达量差异分析部分，我们使用edgeR软件包进行分析。extractExpr和DGE是独立编写的分析脚本。 命令行： cd ~/RNA-seq/ extractExprdata/genes.gtf C460_expression/isoforms.fpkm_tracking 运行完后，在当前目录[RNA-seq目录]下，生成文件夹C460-vs-P460。 P460_expression/isoforms.fpkm_tracking DGE -d expression.xls -g1 C460 -g2 P460 -o differential 可看到生成如下结果：

参数说明：-p 1 表示使用一个线程; -o 表示生成的文件夹名称; genes.gtf是cuffdiff进行差异分析的必需基因结构文件。

结果文件说明：分别产生了cds, gene, isoform,tss的四套差异表达结果，.diff表示差异结果的文件；fpkm_tracking表示表达量估算的文件；count_tracking表示表达量计数的文件。

6. 使用DAVID数据库进行功能搜索

DAVID是一个免费数据库，并提供基因功能搜索，分析和展示的功能，可以分析多种格式，不同物种的基因数据（http://david.abcc.ncifcrf.gov）。按照网址访问DAVID网站，在Start Analysis标签下输入基因编号的列表，选择对应的编号格式，提交即可完成相关搜索和分析。

四、 补充：软件安装和使用环境设置

在计算过程中，常需要自己安装软件。在这补充部分中，主要是关于bowtie,tophat,cufflinks软件的安装和使用环境设置。

首先，分别在bowtie官网[http://bowtie-bio.sourceforge.net/bowtie2/index.shtml]；tophat官网[http://tophat.cbcb.umd.edu/]；cufflinks官网[http://cufflinks.cbcb.umd.edu/]下载对应的软件，例如：unzip bowtie2-2.1.0-linux-x86_64.zip；tophat-2.0.10.Linux_x86_64.tar.gz；cufflinks-2.1.1.Linux_x86_64.tar.gz等软件。

之后，上传软件到服务器，解压文件，把路径导入到环境中，如在RNA-seq目录下，新建bin文件夹，把软件上传到bin文件夹下：

命令行：

cd RNA-seq||mkdirRNA-seq;cd RNA-seq;mkdir bin;cd bin #在RNA-seq下建立文件夹bin #使用filezilla把软件上传到RNA-seq/bin目录下 unzip bowtie2-2.1.0-linux-x86_64.zip export PATH=$HOME/RNA-seq/bin/bowtie2-2.1.0:$PATH tarxvzf tophat-2.0.10.Linux_x86_64.tar.gz export PATH=$HOME/RNA-seq/bin/tophat-2.0.10.Linux_x86_64:$PATH tarxvzf cufflinks-2.1.1.Linux_x86_64.tar.gz export PATH=$HOME/RNA-seq/bin/cufflinks-2.1.1.Linux_x86_64:$PATH 版本核查：

可以看到，环境中的bowtie2版本是2.1.0, tophat版本是2.0.10, cufflinks是2.1.1，正是我们安装的软件。

本文来源：https://www.bwwdw.com/article/g4q5.html