RNA-seq流程

RNA-seq能够更加详细地刻画不同病理或生理状态下转录组的改变。基于RNA-seq数据在不同状态间进行差异表达基因的识别是研究疾病机制以及临床应用的主要手段。
流程参考了生物信息课堂上老师给的文档。
数据参考文章https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50177
数据选了两个进行试验:SRR957677 SRR957678
软件篇:
1.Aspera:一款高速传输软件,可以帮助我们快速下载序列文件(安装教程参考这篇文章https://www.jianshu.com/p/19dcff88a1eb

 wget https://download.asperasoft.com/download/sw/connect/3.6.2/aspera-connect-3.6.2.117442-linux-64.tar.gz
 tar -zvxf aspera-connect-3.6.2.117442-linux-64.tar.gz
 sh aspera-connect-3.6.2.117442-linux-64.sh
 ~/.aspera/connect/bin/ascp -h
4.JPG

2.安装R语言,在R语言中安装DESeq2,用DESeq2进行基因表达差异分析(参考文章https://www.jianshu.com/p/4d0812195b65

 sudo apt install r-base-core
8.JPG
R
#进入R语言包
> source("https://bioconductor.org/biocLite.R")
>biocLite("DESeq2")

3.hisat2:比对基因组工具

wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip
unzip hisat2-2.1.0-Linux_x86_64.zip
echo 'export PATH=~/hisat2:$PATH'>>~/.bashrc
 source ~/.bashrc

7.JPG

4.BWA:序列比对软件

wget https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.15.tar.bz2
tar -jxvf bwa-0.7.15.tar.bz2
cd bwa-0.7.15
make

1.JPG

5.Samtools:samtools是一个用于操作sam和bam文件的工具合集(安装参考https://www.jianshu.com/p/0e05b471189f)。
6.htseq-count利用 htsep-count 计算比对到每个基因的短序列数目(安装参考https://www.jianshu.com/p/5750e8e6fd7e
7.fastp: 一款超快速全功能的FASTQ文件自动化质控,过滤,校正,预处理软件,可参考fastp说明文件 https://github.com/OpenGene/fastp
conda install -c bioconda fastp
5.JPG
6.JPG

1.数据获取
利用Aspera下载原始数据时(语句如下wget anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/srainstant/reads/ByRun/sra/SRR/SRR957/SRR957677/SRR957677.sra)出现下图问题:
2.JPG

加参数 -P 33011还是同样的问题,利用prefetch SRR957677同样显示via http....,最后只好在ENA上找到序列的fastq文件,利用wget下载,唯一问题就是速度很慢特别耗时间。

wget https://www.ncbi.nlm.nih.gov//sra/sra-instant/reads/ByRun/sra/SRR/SRR957/SRR957677/SRR957677.sra
wget https://www.ncbi.nlm.nih.gov//sra/sra-instant/reads/ByRun/sra/SRR/SRR957/SRR957678/SRR957678.sra

速度太慢了,一秒十几kb,最后复制ENA上的网址,用迅雷下了,也不是很快一秒一二百kb,总比wget好一点(中间还试了wget anonftp@ftpprivate.ncbi.nlm.nih.gov:/sra/srainstant/reads/ByRun/sra/SRR/SRR957/SRR957677/SRR957677.sra也没成功)。

3.JPG

下载完直接拖进MobaXterm文件列表就行。

fastp -i /home/aa/2/SRR957677.fastq.gz -I /home/aa/2/SRR957678.fastq.gz -o SRR957677.clean.fastq -O SRR957678.clean.fastq --trim_front1=5 --trim_front2=5 --trim_tail1=1 --trim_tail2=1  --html -w 1

生成去除低质量序列的.clean.fastq 同时生成可视化文件 fastp.html 可查看各项质量指标。


4.JPG

3.去除重复性的冗余序列

remove_redunt_pair_fsq.py SRR957677.clean.fastq SRR957678.clean.fastq 1 1unique SRR957677.clean.uniq.fastq SRR957678.clean.uniq.fastq

生成.clean.uniq.fast格式的文件。
4.利用比对软件 bwa 将短序列比对到参考基因组

bwa mem -t 2 -M /home/aa/2/hg19.exons.gtf SRR957677.clean.fastq > SRR957677.sam
bwa mem -t 2 -M /home/aa/2/hg19.exons.gtf SRR957678.clean.fastq > SRR957678.sam

生成sam格式文件。
5.利用 samtools 对 sam 格式的比对文件进行处理,以便进行后续分析。

samtools view -bt /home/aa/2/hg19/hg19.fa.fai -@ 2 -o SRR957677.bam SRR957677.sam 2>>samtools.log
samtools view -bt /home/aa/2/hg19/hg19.fa.fai -@ 2 -o SRR957678.bam SRR957678.sam 2>>samtools.log

为排序后的bam文件建立索引

samtools index SRR957677.sort.bam
samtools index SRR957678.sort.bam

利用samtools对排序后的bam文件进行可视化

samtools tview SRR957677.sort.bam /home/aa/2/hg19/hg19.fa
samtools tview SRR957678.sort.bam /home/aa/2/hg19/hg19.fa

6.利用 htsep-count 计算比对到每个基因的短序列数目

htseq-count -f sam -s no -t CDS -i ID -m union --nonunique=none --secondary-alignments=ignore -o SRR957677_assigned.sam SRR957677.sam /home/aa/hg19/hg19.gff > SRR957677_expression.counts
htseq-count -f sam -s no -t CDS -i ID -m union --nonunique=none --secondary-alignments=ignore -o SRR957678_assigned.sam SRR957678.sam /home/aa/hg19/hg19.gff > SRR957678_expression.counts

生成文件:SRR957677_expression.counts SRR957678_expression.counts

记录比对到每个基因上的短片段数目, 代表基因表达水平的高低。

©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容