10月20日 RNA-seq 实战(一)服务器没网 用户没root权限

简单地记录一下我的一次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语言来做可视化

后记:

这个流程现在看看是很简单,但是从一点也不会到跑下来遇到的问题都是很困难的,特别是没有人带你,还是要感谢各路大神,各种论坛。 说不定这个流程在我服务器上跑是可以的,但是别人的就不行,所以多多交流哈

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

推荐阅读更多精彩内容