搜索
查看: 901|回复: 14

[mRNA-seq] 跟随群主博客学转录组分析

[复制链接]

4

主题

8

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
QQ
发表于 2017-1-3 11:41:38 | 显示全部楼层 |阅读模式
拖延症啊拖延症,,, 谁来救救我,,救命啊,,,,

根据群主Jimmy的博客内容进行的一个转录组分析的过程化复制,原文链接:http://www.bio-info-trainee.com/2218.html

首先下载相关内容,请移步群主博客,我不做搬运工。
基本流程如下:

fastqc
hisat2 比对
统计比对信息
表达量定量——htseq-count
制作基因表达列表
用whatever进行差异表达分析。

质控数据群主做了,表示没有异常,我,,,懒。
hisat2比对按照manual进行, 先index,然后比对, 记得输出一下log,log里有记录map信息。
示例代码(抄袭群主):~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR957677.fastq -S control_1.sam 2>control_1.log
这里群主是先将基因组(人的hg19.fa)所在的绝对路径保存到了一个变量$reference里,很巧妙的一个替代。 我这里是用ln -s 做了一个软连接,相当于快捷方式,一个道理。可以将这四个代码写成一个循环,或者写到一个.sh脚本里丢给后台:


比对得到的是sam文件,之后要使用samtools里的功能来进行排序,统计跟格式转换等。
在此可能由于版本问题,博客里的命令行pipe我就是用不了,所以还是自己写shell 慢慢来的。
samtools view -bS XXX.sam > XXX.bam
samtools sort -n -@ 4 XXX.bam XXX.sorted

排序完了就可以进行统计了,还是用不了博客里的pipe命令,略郁闷,自己写shell,或者一个一个来,反正就4个。
htseq-count -f bam -s no -i gene_name control_1.sorted.bam gencode.v25lift37.annotation.gtf >control_1.genecounts 2>control_1.htseq.log

得到四个genecounts文件后,就可以制作表达列表了。
笨方法: 先用awk将每个文件里面的counts这一列提出来,然后再用paste命令合并。

最后用edgeR进行差异分析,使用edgeR进行差异表达分析-前期摸索版
http://www.biotrainee.com/forum. ... 607&fromuid=728
(出处: 生信技能树)

最后的结果跟原文章比对一下就好。



上一篇:致病性SNP分布在基因什么位置?
下一篇:Perl-000【作业示例】
Yes, I am serious.
回复

使用道具 举报

4

主题

21

帖子

285

积分

版主

Rank: 7Rank: 7Rank: 7

积分
285
发表于 2017-1-3 13:25:40 | 显示全部楼层
厉害了我的哥
回复 支持 反对

使用道具 举报

1

主题

15

帖子

116

积分

注册会员

Rank: 2

积分
116
发表于 2017-1-10 15:01:37 | 显示全部楼层
为什么不用stringtie+ballgown分析?
回复 支持 反对

使用道具 举报

0

主题

6

帖子

60

积分

注册会员

Rank: 2

积分
60
发表于 2017-1-12 10:58:22 | 显示全部楼层
请教,HISAT比对转录组数据不需要用gtf文件来指导吗
回复 支持 反对

使用道具 举报

515

主题

907

帖子

3067

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3067
发表于 2017-1-13 09:38:38 | 显示全部楼层
lipc 发表于 2017-1-12 10:58
请教,HISAT比对转录组数据不需要用gtf文件来指导吗

比对的时候,不需要,counts的时候才需要
回复 支持 反对

使用道具 举报

515

主题

907

帖子

3067

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3067
发表于 2017-1-13 09:39:01 | 显示全部楼层
starhqy2 发表于 2017-1-10 15:01
为什么不用stringtie+ballgown分析?

人类的转录组,不需要搞这些东西,简单的看看表达量就好了
回复 支持 反对

使用道具 举报

515

主题

907

帖子

3067

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3067
发表于 2017-1-13 09:39:20 | 显示全部楼层

你的哥的哥在这里
回复 支持 反对

使用道具 举报

0

主题

6

帖子

60

积分

注册会员

Rank: 2

积分
60
发表于 2017-1-13 22:58:34 | 显示全部楼层
本帖最后由 lipc 于 2017-1-13 22:59 编辑
Jimmy 发表于 2017-1-13 09:38
比对的时候,不需要,counts的时候才需要

请教一下,用tophat2比对转录组时,构建完索引之后,把基因组数据删掉后,比对提示如下,比对可以正常进行
[2017-01-13 22:31:39] Checking for Bowtie index files (genome)..
[2017-01-13 22:31:39] Checking for reference FASTA file
        Warning: Could not find FASTA file

保留基因组数据再比对一次。
两次比对的align_summary.txt,信息好像是一致的。 第一次的提示是不是不用管。构建好索引后,没有fasta的基因组文件是不是不影响结果?





回复 支持 反对

使用道具 举报

515

主题

907

帖子

3067

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3067
发表于 2017-1-14 21:14:22 | 显示全部楼层
lipc 发表于 2017-1-13 22:58
请教一下,用tophat2比对转录组时,构建完索引之后,把基因组数据删掉后,比对提示如下,比对可以正常进行 ...

tophat2我不用的,推荐换hisat2
回复 支持 反对

使用道具 举报

0

主题

6

帖子

60

积分

注册会员

Rank: 2

积分
60
发表于 2017-1-15 16:24:15 | 显示全部楼层
Jimmy 发表于 2017-1-14 21:14
tophat2我不用的,推荐换hisat2

开始先试了一下tophat2,个人电脑跑不动,现在用hisat2,但是群主博客里对SAM文件操作的代码确实用不了,不知道什么原因
[Bash shell] 纯文本查看 复制代码
ls *sam |while read id;do (nohup samtools sort -n -@ 5 -o ${id%%.*}.Nsort.bam $id  );done
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|关于我们|手机版|小黑屋|生信技能树    

GMT+8, 2017-4-26 00:42 , Processed in 0.034269 second(s), 33 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.