简单地记录一下我的一次RNA-seq实战,由于我们实验室的服务器没有联网,我之前一直在自己的PC上,安装虚拟机或者用OS系统跑数据,出现许多问题,而且速度很慢,后来询问了洲更大神,直接按照biostar handbook 来设置我的电脑Mac,但是我又遇到了很多问题,不过biostar handbook 确实很不错,我后续的学习应该要按照这本书来,但是这次时间比较急,要赶紧先跑一遍,学会大致流程再说。那我就跌跌撞撞根据百度,文献,问大神的模式开始学习了。
整个学习过程是完全没有逻辑和顺序的,但是我现在是整理好的顺序,这样后面的人来看操作起来就很方便了。
1安装软件的过程
先要安装这些软件,这些软件直接在百度搜索,然后去官网下载软件包,Linux 下直接解压可以用的就可以,建议不要源码编译,我源码编译了一下,但是服务器没有联网,缺少lib 库 反正没有直接下载软件包来的方便。
fastqc
fastx_toolkit
trimmomatic
bowtie2
tophat
cufflinks
1传输文件ssh
scp -P 22 -r 要传输的文件路径 will@10.10.10.10 你要保存的文件路径
2解压.gz .tar.gz
解压这些文件,解压不同类型的方法不一样,具体可以参考这篇文章:
http://www.cnblogs.com/qq78292959/archive/2011/07/06/2099427.html
*.tar 用 tar -xvf 解压
*.gz 用 gzip -d或者gunzip 解压
*.tar.gz和*.tgz 用 tar -xzf 解压
*.bz2 用 bzip2 -d或者用bunzip2 解压
3 把文件PATH修改
1 cd 回到主目录
2 ls -all 查看所有的文件 找到 .bash_profile ( 可能不同系统命名不一样)
3 vim .bash_profile
增加 exportPATH=$PATH:/public/home/你的用户名/rna-seq/FastQC/
4 FASTqc使用指南检测测序的质量
这个看我之前写的一篇文章专门记录fastqc 以及遇到的问题 或者参考这篇https://zhuanlan.zhihu.com/p/20731723
主要命令如下:
1 主文件下 chmod 755 fastqc
2 fastqc -o result -t 8 gene_1.fq
然后会生成一个报告文件,具体分析可以看之前的文章
5 Trimmomatic处理reads 可以参考
//www.greatytc.com/p/36891a89ed6e
http://blog.sina.com.cn/s/blog_4b91a9e50101ock4.html
主要就是java 的程序
java -jar trimmomatic-0.36.jar PE -threads20 -phred33 -trimlog logfile reads1.fq reads2.fq reads1.clean.fq reads1.unpaired.fq reads2.clean.fq reads2.unpaired.fq ILLUMINACLIP:/path/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10\ LEADING:3TRAILING:3SLIDINGWINDOW:4:15MINLEN:50
PE/SE 设定对Paired-End或Single-End的reads进行处理,其输入和输出参数稍有不一样。 -threads 设置多线程运行数 -phred33 设置碱基的质量格式,可选pred64ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 切除adapter序列。参数后面分别接adapter序列的fasta文件:允许的最大mismatch 数:palindrome模式下匹配碱基数阈值:simple模式下的匹配碱基数阈值。 LEADING:3 切除首端碱基质量小于3的碱基 TRAILING:3 切除尾端碱基质量小于3的碱基 SLIDINGWINDOW:4:15 Perform a sliding window trimming。Windows的size是4个碱基,其平均碱基 质量小于15,则切除。 MINLEN:50 最小的reads长度 CROP: 保留reads到指定的长度 HEADCROP: 在reads的首端切除指定的长度 TOPHRED33 将碱基质量转换为pred33格式 TOPHRED64 将碱基质量转换为pred64格式
会出现4个文件
http://blog.sciencenet.cn/blog-2458445-933838.html
paired 的一对 unpaired 的一对
用biotie2 处理paired 的一对
其实我还是有点不理解,下次理解了再补上这里的
6 bowtie2 (map的常用工具 )(短序列拼接至模版基因组)
好了现在reads处理好了,接下来要map了
1 先要建立参考基因组的索引 对参考基因组fasta 格式文件建立index ( 去网上下载参考基因组的fa文件)
bowtie2-build genome.fa genome
然后你会看到 genome.1.bt2 4个 还有 genome.rev.1.bt2 2 个
7 tophat通过调用bowtie2来完成map
tophat -p 30 -G annotation.gff3 -o qc-thout genome qc-1.fq qc-2.fq
( annotation.gff3 为基因的注释信息 genome为基因组的序列,qc-1.fq qc-2.fq为质控后的双端测序数据)
出现错误一 (注意啦 注意啦)
Error: Could not find Bowtie 2 index files (genome.*.bt2)
这个是网上找的,和我基本一样,我也是按照网上攻略输入的语句,但是就是出现问题。然后我捣鼓了好久,去了https://www.biostars.org/p/203950/这个biostar 网站,还是不错的,不过里面的英文有些理解不了,要参考一下中文的. 解决如下:
1 最好把要跑的文件和参考基因组的文件都放在一个文件夹下面
2 (我的问题主要是这个) basename 很重要, 一定要把basename 改成和基因的注释文件一样的开头,不能下载下来是什么就用什么
3 我把基因的注释文件后面的 .fa拿掉了,发现才可以跑起来,不知道为什么
出现错误二 ( 注意啦注意啦)
跑了大约1个小时的时候突然断了
Tophat Error "segment-based junction search failed"
这是咋搞得呢,我查了查https://www.biostars.org/p/146557/
这篇英文说的还是挺清楚的 我线程用的太多了 我用了 80.....
好吧我把 -P 调节到30 就解决了
OK 跑完出现了几个文件:
accepted_hits.bam align_summary.txt deletions.bed insertions.bed junctions.bed logs prep_reads.info unmapped.bam
参考一下这个解析很清楚:http://www.360doc.com/content/17/0802/19/45954458_676165658.shtml
8 cufflinks查看样本间基因的差异表达计算
主要用于基因表达量的计算和差异表达基因的寻找。
Cufflinks程序主要根据Tophat的比对结果,依托或不依托于参考基因组的GTF注释文件,计算出(各个gene的)isoform的FPKM值,并给出trascripts.gtf注释结果(组装出转录组)。
1 用 cufflinks 对Tophat比对的每个结果(bam文件)进行组装
cufflinks -p 30 -u -g annotation.gff3 --library-type fr-unstranded - o S3 accepted_hits.bam
出现下面几个文件
genes.fpkm_tracking isoforms.fpkm_tracking skipped.gtf transcripts.gtf
2 将 Cufflinks 组装出来的转录本注释文件(gtf)结尾的 路径汇总到文件 assemblies.txt
然后用cuffmerge将各个样本的组装结果与下载的基因注释文件整合在一起,创建一个整体的注释信息:
cuffmerge -p 30 -o merge_out -g annotation.gff3 -s genome.fa assemblies.txt
注意这个assemblies.txt文件要自己建立,然后写入路径
3 最后啦啦啦
用cuffdiff -o diff-out -b genome.fa -p 40 -L S1,G1 -u merged.gtf accepted_hits.bam
9未完待续后面用R语言来做可视化
后记:
这个流程现在看看是很简单,但是从一点也不会到跑下来遇到的问题都是很困难的,特别是没有人带你,还是要感谢各路大神,各种论坛。 说不定这个流程在我服务器上跑是可以的,但是别人的就不行,所以多多交流哈