如何利用linux分析转录组数据库,对转录组测序数据进行分析以及注释

随着二代测序技术的高速发展，人们获得了大量的转录组数据序列如何从数据中挖掘具有生物意义的信息已经成为很多研究的关键，对未知基因的功能进行预测和注释就是一个重要问题

这篇文章主要是跟着 刘粉香，杨文国，孙勤红，三位老师的文献对测序数据进行分析以及GO注释，旨在学习测序数据分析的方法以及GO注释的方法

数据来源于NCBI上的SRA数据库 SRR063784

1.原始测序数据下载与质控

#下载数据

~/.aspera/connect/bin/ascp -T -i /home/qiujunhui/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -l 200m anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR063/SRR063784/SRR063784.sra ./

#解压

fastq-dump SRR063784.sra

#质控

mkdir fastqc_out

fastqc -o fastqc_out/ --noextract SRR063784.fastq

#用Trimmomatic切除reads的phred值小于20的

java -jar ~/Biosofts/Trimmomatic-0.38/trimmomatic-0.38.jar SE -phred33 SRR063784.fastq SRR063784_clean.fastq LEADING:20 TRAILING:20

TrimmomaticSE: Started with arguments:

-phred33 SRR063784.fastq SRR063784_clean.fastq LEADING:20 TRAILING:20

Automatically using 4 threads

Input Reads: 31028774 Surviving: 26555513 (85.58%) Dropped: 4473261 (14.42%)

TrimmomaticSE: Completed successfully

#可以看到有4473261条reads因为碱基质量过低被舍弃，不做去除接头是因为原始数据没有接头

2.从头拼接

#由于服务器内存太小！！！所以我直接选用最优Kmer值37，原文中是准备了5个

velveth chickpea_37 37 -short -fastq SRR063784_clean.fastq #velveth准备数据

velvetg chickpea_37 -read_trkg yes #velvetg做转录组拼接

#然后再chickpea_37文件夹里面有8个文件，contigs.fa为拼接后的contigs

最后运行oases对Velvet拼接到的contig进行进一步的拼接获得转录本

oases安装与使用见：https://www.plob.org/article/918.html

./oases ~/chickpea_37

运行结束后文件夹chickpea_37中产生两个文件，分别是transcripts.fa和contig-ordering.txt。前者为包含组装得到的transcripts文件，后者记录了每一个transcripts中的contigs组成情况

3.基因注释与功能分类

做到这里我已经连续做了6.5个小时了，可能是因为第一次做吧，对流程不熟悉，软件的安装与使用都要自己摸索。。。头很痛

头痛不是因为做的时间长，而是实验条件有限制，做不了本地blast对转录组数据进行注释，但是放弃又觉得可惜。。。我就接着把后面的教程写下来，等以后有条件再做后续实验

到这一步是要通过比对已知数据库中的同源基因的信息推断未知基因的功能

这里用的软件是NCBI上的blast工具，可以在NCBI上下载，解压

要进行本地Blast搜索，首先需要从NCBI的ftp站点下载并且格式化数据库nr.gz,额大概100G左右(这还是压缩后的大小)，所以就。。。。

1.将下载的nr.gz放在目录ncbi-blast-2.2.25+/bin/中，并解压。利用文件夹bin/中的makeblastdb格式化数据库，代码如下

make-blastdb -in nr -dbtype prot -parse_seqids -out nr.db

in (nr)是待格式化的文件

dbype (prot)给出数据库的类型(蛋白质数据库)

-parse_seqids 启动序列ID解析

格式化数据库后，就可以运行Blast将拼接得到的transcipts比对到本地nr数据库 代码如下

blastx -query transcripts.fa -out transcripts.xml -db ~/Biosofts/blast+/bin/nrdb -outfmt 5 -evalue 1.OE -6 -max_target_seqs 10 -num_threads 20

-query 给出输入带比对的数据文件

-out 指定输出文件

-db 指定本地数据库

-outfmt 指定输入数据格式(xml格式)

-evalue 设置输出结果的E-value

-num_threds 使用多线程运算

4.利用Blast2Go进行GO注释

Blast+只是一种预测新基因功能的基本工具，仅通过Blast无法得到新基因的GO注释信息 可以将Blast的输出文件(xml文件)作为Blast2God的输入文件对其进行GO注释。

然而依然需要配置本地数据库(100G！！！)

1.进入Blast2Go主页(http://www.blast2go.com/)下载并安装软件，在Linux下打开Blast2Go运行界面输入代码

Javaws -Xnosplash blast2go*.jnlp

2.选择合适的数据库后，点击start便开始Blast比对

(1)直接输入xml文件

(2)输入序列文件

(3)选择Blast运行方式

(4)Blast设置

3.进行Mapping

4..进行Annotating

这些都是软件的使用方法，慢慢摸索总会的~