搜索
查看: 1939|回复: 6

[lncRNA-seq] 我的初尝试找lncRNA 希望大家能多多指教

[复制链接]

2

主题

9

帖子

147

积分

注册会员

Rank: 2

积分
147
发表于 2018-2-12 17:58:43 | 显示全部楼层 |阅读模式

  • title: lncRNA
  • author: qliu
  • date: 2018年2月12日
  • output: lncRNA.html

  • 使用Tophat2比对到水稻参考基因组上
tophat2 \--min-intron-length 8 \--max-intron-length 60000 \-p 8 \-r -25 \--library-type fr-firststrand \-o mut/ \--transcriptome-index /public/home/qliu/temp/zc_data/lncRNA/trancript_index/all \/public/home/qliu/zhouchao/bowtie2_index/rice \/public/home/qliu/zhouchao/cleandata/mut_R1.fq \/public/home/qliu/zhouchao/cleandata/mut_R2.fqtophat2 \--min-intron-length 8 \--max-intron-length 60000 \-p 8 \-r -25 \--library-type fr-firststrand \-o WT/ \--transcriptome-index /public/home/qliu/temp/zc_data/lncRNA/trancript_index/all \/public/home/qliu/zhouchao/bowtie2_index/rice \/public/home/qliu/zhouchao/cleandata/WT_R1.fq \/public/home/qliu/zhouchao/cleandata/WT_R2.fq
  • 使用Cufflinks 2.0进行组装
cufflinks -p 8 --library-type fr-firststrand -o mut -g /public/home/qliu/temp/zc_data/lncRNA/trancript_index/all.gff /public/home/qliu/temp/zc_data/lncRNA/tophat2/mut/accepted_hits.bamcufflinks -p 8 --library-type fr-firststrand -o WT -g /public/home/qliu/temp/zc_data/lncRNA/trancript_index/all.gff /public/home/qliu/temp/zc_data/lncRNA/tophat2/WT/accepted_hits.bam
  • 使用Cuffmerge将所有的转录本进行组合(merge),形成一个新的转录组
cuffmerge -p 8 -g /public/home/qliu/temp/zc_data/lncRNA/trancript_index/all.gff  -s /public/home/qliu/temp/zc_data/lncRNA/genome/all.fa  assemblies.txt
  • 使用cuffcompare找新转录本
## 当class_code为”u”时,为lincRNA,当class_code为”i”时,为intronic lncRNA,当class_code为”x”时,为anti-snselncRNA。cuffcompare -o merged_lncRNA -r /public/home/qliu/temp/zc_data/lncRNA/trancript_index/all.gff -p 4 /public/home/qliu/temp/zc_data/lncRNA/cuffmerge/merged_asm/merged.gtf
  • 使用Cuffdiff基于新形成的转录组去估计所有转录本的丰度
cuffdiff -p 4 -o diffOut \--library-type fr-firststrand \-b /public/home/qliu/temp/zc_data/lncRNA/genome/all.fa \-L WT,mut \-u \/public/home/qliu/temp/zc_data/lncRNA/cuffmerge/merged_asm/merged.gtf \/public/home/qliu/temp/zc_data/lncRNA/tophat2/WT/accepted_hits.bam \/public/home/qliu/temp/zc_data/lncRNA/tophat2/mut/accepted_hits.bam文献操作流程
  • Tophat2比对形成的BAM文件
  • 去除掉所有没有链信息的转录本和在其他转录本正链上所有在500bp范围内的单个外显子的转录本
  • 将覆盖在已知的mRNAs(包括TE相关的基因)和转录本FPKM值小于0.5单个外显子转录本为2)并且转录本长度小于200bp的转录本去除掉
  • 过滤后得到的转录本即为新的长的转录本
  • 然后使用CPC软件预测可能能编码的转录本
  • 所有CPC打分大于0的都被去除
  • 剩下的转录本使用HMMER软件分析提取在Pfam数据库中包含已知蛋白结构域的转录本
  • 得到的转录本即为可靠的表达的lncRNA
  • 使用CummeRbund 中的R包 csCpecificity() 计算每个转录本的组织特异性打分
  • 结直肠癌lncRNA pipeline
TopHat/Cufflinks/CummeRbund使用介绍
awk '$22 ~ /j/ { print }' cuffcompare_combined.gtf > Alternatively_spliced.gtfawk '{ if ($5-$4>200) print $0 }'  Alternatively_spliced.gtf > Alternatively_spliced_200.gtf##采用excel去除exon_numbers <= 1 得到 Alternatively_spliced_200_exon2.gtf## bedtools getfasta [OPTIONS] -fi <fasta> -bed <bed/gff/vcf> -fo <fasta> -s## 注意-s参数bedtools getfasta -fi /public/home/qliu/temp/zc_data/lncRNA/genome/all.fa -bed Alternatively_spliced_200_exon2.gtf -fo lncRNA.fa -sCPC在线预测 : 剔除值CODING POTENTIAL SCORE < 0的
CPC
提交上面得到的lncRNA.fa文件
task ID : 903784E0-0FF2-11E8-9838-B4256C0D250A
Pfam蛋白编码能力预测所需软件安装hmmpress Pfam-A.hmm## 这一步需要剔除重复的一些坐标,由于存在重复## 注意这里的lncRNA.fa应该是CPC预测后的fa文件,但是我这里没有剔除。./pfam_scan.pl -fasta /public/home/qliu/temp/zc_data/lncRNA/cuffcompare/lncRNA/lncRNA.fa -dir /public/home/qliu/bio_soft/pfam_data/ -as -outfile /public/home/qliu/temp/zc_data/lncRNA/cuffcompare/lncRNA/result.fa#参数说明:-dir Pfam_data_file_dir 包含Pfam数据文件的目录[必须]-fasta fasta_file 包含序列的输入文件名 [必须]-outfile output_file 输出文件名 [不指定则输出在命令行中]-as 预测Pfam-A数据库匹配的activesites[默认关闭]
pfamscan 在线预测检测
Pfamscan

致谢所有能把自己的学习经验分享出来的小伙伴,才能让更多向我这样的小白能更快的入门。




上一篇:问一个大家bam文件修改HI标签的问题
下一篇:Linux服务器非root安装mysql
回复

使用道具 举报

2

主题

9

帖子

147

积分

注册会员

Rank: 2

积分
147
 楼主| 发表于 2018-2-12 18:11:17 | 显示全部楼层
文献链接lncRNA_2014_rice_GB
回复 支持 反对

使用道具 举报

23

主题

37

帖子

344

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
344
发表于 2018-2-14 19:51:07 | 显示全部楼层
qliu 发表于 2018-2-12 18:11
文献链接lncRNA_2014_rice_GB

写的很棒,插入代码,论坛首页有教程,会看起来比较好一些。
回复 支持 反对

使用道具 举报

2

主题

9

帖子

147

积分

注册会员

Rank: 2

积分
147
 楼主| 发表于 2018-2-15 00:13:13 | 显示全部楼层
张翼翔 发表于 2018-2-14 19:51
写的很棒,插入代码,论坛首页有教程,会看起来比较好一些。

本来就是markdown转过来的 不知道为什么没成功。。
回复 支持 反对

使用道具 举报

0

主题

3

帖子

61

积分

注册会员

Rank: 2

积分
61
发表于 2018-4-14 21:47:47 | 显示全部楼层
个人认为你的分析有两个问题
回复 支持 反对

使用道具 举报

0

主题

3

帖子

61

积分

注册会员

Rank: 2

积分
61
发表于 2018-4-14 21:51:35 | 显示全部楼层
1)bedtools getfasta -fi /public/home/qliu/temp/zc_data/lncRNA/genome/all.fa -bed Alternatively_spliced_200_exon2.gtf -fo lncRNA.fa

这个只能拿到从start 到end 间的所有序列 不能去除内含子序列
2)awk '$22 ~ /j/ { print }' cuffcompare_combined.gtf > Alternatively_spliced.gtfawk '{ if ($5-$4>200) print $0 }'  Alternatively_spliced.gtf > Alternatively_spliced_200.gtf#
这个也有点问题$5-$4>200 不对 也是终止减去起始的 没有考虑中间内含子长度,仔细看看gtf文件
回复 支持 反对

使用道具 举报

2

主题

9

帖子

147

积分

注册会员

Rank: 2

积分
147
 楼主| 发表于 2018-5-9 10:28:58 | 显示全部楼层
charleslywang 发表于 2018-4-14 21:51
1)bedtools getfasta -fi /public/home/qliu/temp/zc_data/lncRNA/genome/all.fa -bed Alternatively_spli ...

这个还真没考虑欸。。谢谢谢谢了
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2018-10-24 01:19 , Processed in 0.111971 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.