干货 | 一文教会你如何采用Linux系统处理RNAseq测序数据

    转录组是指细胞在某一功能状态下转录出来的所有RNA的总和。转录组测序(Transcriptome sequencing)是基于Illumina HiSeq测序平台检测细胞内所有mRNA的一项技术,能够快速获得细胞在某一状态下所有的转录本信息,因而被广泛应用于基础研究、药物研发和临床诊断等多个领域。

    在获得原始测序数据(FASTQ文件格式)后,该如何就FASTQ测序文件进行后续分析处理呢?下面,让我们一起来学习吧!

前期准备:

硬件:

Linux系统,>4G运行内存

软件:

Miniconda、FastQC、Cutadapt、Hisat2、SAMtools、subread(FeatureCounts)

PART1 测序数据质量评估采用FastQC软件进行分析

    FastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)是一款基于Java的软件。无论是Windows/Linux系统运行FastQC,首先都必须安装java(注意Java的版本不能超过1.6,否则FastQC会无法运行)。其作用是对测序数据进行质量评估。

FastQC安装:

$ cd path/to/fastqc

$ chmod 755 fastqc

$ ./fastqc

$ sudo ln -s /path/to/FastQC/fastqc /usr/local/bin/fastqc

运行格式:

$ fastqc -o outdir -t threads fastq1 fastq2 ...

主要参数:

-o --outdir FastQC 所生成的报告文件的储存路径

-t --threads 程序运行的线程数

示例:

$ fastqc -o FASTQC/ -t 8 Ctrl-1_combined_R1.fastq.gz Ctrl-1_combined_R2.fastq.gz Ctrl-2_combined_R1.fastq.gz Ctrl-2_combined_R2.fastq.gz Ctrl-3_combined_R1.fastq.gz Ctrl-3_combined_R2.fastq.gz Treated-1_combined_R1.fastq.gz Treated-1_combined_R2.fastq.gz Treated-2_combined_R1.fastq.gz Treated-2_combined_R2.fastq.gz Treated-3_combined_R1.fastq.gz Treated-3_combined_R2.fastq.gz

$ multiqc ./

运行结果:

运行结束后,生成两个文件:.html网页文件及.zip压缩文件。

PART2 测序数据过滤采用Cutadapt软件进行处理

    测序过程中少量reads会测到接头序列,或者测序长度过长时活导致3’端碱基质量过低,因此需要对原始数据进行预处理。采用Cutadapt(https://cutadapt.readthedocs.io/en/stable/installation.html)可以帮助我们:(1)去除接头序列;(2)去除5’或3’末端质量值较低或含N的碱基;(3)去除平均质量值低于30的序列;(4)去除trim后reads长度过短的序列。

Cutadapt安装(需提前安装miniconda):

$ conda install -c bioconda cutadapt

运行格式:

$ cutadapt -a ADAPTER_FWD -A ADAPTER_REV -o out.1.fq -p out.2.fq reads.1.fq reads.2.fq

主要参数:

-a 正向接头序列(单端测序时仅有该参数)

-A 反向接头序列(双端测序时增加该参数)

-q 表示最低质量值,将低于此数值的碱基去除,一般设为30

-m 去除去接头后短于此数值的reads

--trim-n 去除含N的碱基

-o reads.1.fq去接头后的输出文件

-p reads.2.fq去接头后的输出文件

示例:

$ cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT -q 30 -m 75 --trim-n --report=minimal -o Ctrl-1_out1_R1.fastq.gz -p Ctrl-1_out1_R2.fastq.gz Ctrl-1_combined_R1.fastq.gz Ctrl-1_combined_R2.fastq.gz

运行结果:

获得过滤后的Clean Data。此时可再次运行一次FastQC软件查看过滤后的数据质量。

PART3 采用Hisat2软件与参考基因组进行比对分析

    获得Clean data后,即可与参考基因组进行比对分析,我们采用Hisat2软件(http://ccb.jhu.edu/software/hisat2/index.shtml)进行短reads的比对,以人类参考基因组为例。

Hisat2安装

$ unzip hisat2-2.0.1-beta-OSX_x86_64.zip

$ cp hisat2 hisat2-align-s hisat2-align-l hisat2-build hisat2-build-s hisat2-build-l hisat2-inspect hisat2-inspect-s hisat2-inspect-l $HOME/bin

$ vi ~/.bashrc

$ export PATH=/lustre/home/lcn/chenwen/bin/hisat2-2.0.1-beta:$PATH

$ source ~/.bashrc

构建索引:

i)提取剪接位点及外显子信息

$ extract_splice_sites.py genes.gtf >genome.ss

$ extract_exons.py genes.gtf >genome.exon

ii)创建HISTA2 index

$ hisat2-build -p 20 --s genome.ss --exon genome.exon genome.fa genome_tran

    创建人基因组的索引往往需要很大的内存(200G RAM),因此不建议大家自己创建,常用的物种均可以去下面网站下载现成的(http://ccb.jhu.edu/software/hisat2/index.shtml)。

运行格式:

$ hisat2 -p 8 --dta -x /path/to/file/hg19/genome -1 Ctrl-1_out_R1.fastq.gz -2 Ctrl-1_out_R2.fastq.gz -S Ctrl-1_out.sam

主要参数:

-p 线程数

-x 指定索引文件

-1 指定第一个FASTQ文件

-2 指定第二个FASTQ文件

-S 输出的SAM文件

运行结果:

获得比对后的SAM文件。

PS:若想将测序文件与自己的序列进行比对,可以参考小编之前的推送。

如何将转录组数据mapping到自己的序列并可视化?(HISAT2+Samtools+IGV)

PART4 SAMtools软件将SAM文件转化为BAM文件

    采用SAMtools软件(http://samtools.sourceforge.net)进行sort及格式转化,sort之后文件占位更小,其转化为BAM二进制文件。

SAMtools安装:

$ tar jxvf samtools-0.1.19.tar.bz2

$ cd samtools-0.1.19

$ make

$ cd.

$ cp samtools-0.1.19/samtools $HOME/bin

运行格式:

$ samtools sort -@ 8 -o Ctrl-1_out.bam Ctrl-1_out.sam

运行结果:

获得BAM文件,该结果可采用IGV进行可视化。

PART5 采用FeatureCounts软件进行reads计数

    FeatureCounts是subread软件包中的一个工具,尽管短序列比对工具subread没有bwa和hisat2流行,但其中FeatureCounts工具却应用广泛。其作用是将比对后的结果进行计数,即计算每个区域比对到的reads数。之后对reads数进行归一化处理,计算RPKM,FPKM,TPM等,然后进行差异比较分析。

Subread安装(直接找到官网下载解压即可):

$ wget https://jaist.dl.sourceforge.net/project/subread/subread-1.6.0/subread-1.6.0-Linux-x86_64.tar.gz

$ tar -zxvf subread-1.6.0-Linux-x86_64.tar.gz

运行格式:

$ featureCounts -T 5 -p -t exon -g gene_id -a /path/to/file/genes.gtf -o Ctrl-1.featureCounts.txt /path/to/file/Ctrl-1_out.bam

    其中参考基因组需自行从UCSC下载。

主要参数:

-T 多线程数

-p 针对双端测序数据

-t 指定注释文件中的功能类型,默认为外显子exon

-g 指定注释文件中的属性信息,默认为gene_id

-a 注释文件,GTF或GFF格式文件,区分外显子区域

-o 输出结果文件,同时会输出以该名称结尾的.summary统计文件

运行结果:

得到rawcounts文件。之后可在R中进行后续差异比较分析。

参考文献:

1.    Mihaela Pertea, Daehwan Kim, Geo M Pertea, Jeffrey T Leek, Steven L Salzberg. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Prot 2016. 11(9):1650-1667.

2.    Yang Liao, Gordon K. Smyth, Wei Shi. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2014. 30(7):923-930.

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 211,884评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,347评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,435评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,509评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,611评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,837评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,987评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,730评论 0 267
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,194评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,525评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,664评论 1 340
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,334评论 4 330
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,944评论 3 313
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,764评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,997评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,389评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,554评论 2 349

推荐阅读更多精彩内容