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
2.安装R语言,在R语言中安装DESeq2,用DESeq2进行基因表达差异分析(参考文章https://www.jianshu.com/p/4d0812195b65)
sudo apt install r-base-core
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
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
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
1.数据获取
利用Aspera下载原始数据时(语句如下
wget anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/srainstant/reads/ByRun/sra/SRR/SRR957/SRR957677/SRR957677.sra
)出现下图问题:加参数 -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
也没成功)。
下载完直接拖进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 可查看各项质量指标。
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
记录比对到每个基因上的短片段数目, 代表基因表达水平的高低。