small RNA学习(四):鉴定&注释产生小RNA的基因

1. 使用ShortStack鉴定,注释产生小RNA的基因

1.1 下载安装
$ git clone https://github.com/MikeAxtell/ShortStack.git
$ cd ShortStack/
$ ls
LICENSE  README  ShortStack
1.2 下载依赖的软件(包)

perl, samtools, bowtie & bowtie-build这几个都是平时常用的,不用安装了。
另外的RNALfold & RNAeval(ViennaRNA package)以及einverted (EMBOSS package)以前没接触过,需要安装。见这里

那这几个软件有什么作用呢?
RNALfold & RNAeval:Predicts locally stable RNA secondary structures
einverted:Finds inverted repeats in DNA sequences
bowtie & bowtie-build:建立索引,无gap的短reads比对
samtools:处理sam格式文件

1.3 简单使用
ShortStack \
--readfile /srna_project/data/out.SRR1266859.fastq.gz /srna_project/data/out.SRR1266860.fastq.gz \
--genomefile /srna_project/ref/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa \
--bowtie_cores 6 --sort_mem 10G

会自动生成一个以ShortStack_开头的目录,最终生成如下文件:
$ ls -t | cat
Log.txt
ShortStack_All.gff3
ShortStack_D.gff3
ShortStack_N.gff3
Counts.txt
Unplaced.txt
Results.txt
MIRNAs/
merged_alignments.bam
ErrorLogs.txt

中间会生成out.SRR1266860_readsorted.sam.gz(以queryname排序),out.SRR1266860.bam(以坐标排序)等文件,最后会删除;
参考基因组所在的文件夹下还会生成samtools faidx和bowtie-build产生的索引文件。

2. 结果解读

2.1 ErrorLogs.txt

报错信息(如果有的话);还包含一些运行信息,如bowtie-build, bowtie, samtools sort, and samtools merge的详细运行信息。比如下面的就比较重要:

# reads processed: 20394922
# reads with at least one reported alignment: 16290733 (79.88%)
# reads that failed to align: 1705936 (8.36%)
# reads with alignments suppressed due to -m: 2398253 (11.76%)
Reported 60107848 alignments
# reads processed: 27228269
# reads with at least one reported alignment: 22673414 (83.27%)
# reads that failed to align: 2373241 (8.72%)
# reads with alignments suppressed due to -m: 2181614 (8.01%)
Reported 75813835 alignments
2.2 Results.txt

包含主要结果信息,每一列的具体含义见说明文档。

Locus   Name    Length  Reads   RPM UniqueReads FracTop Strand  MajorRNA    MajorRNAReads   'Complexity DicerCall   MIRNA   PhaseScore' Short   Long    20  21  22  23  24
2.3 Unplaced file

显示了每一条满足--mincov的无法定位的小RNA reads。其中第五列表示能比对到参考基因组上的位置数,>=0的任意数;第一列第二列表示序列和长度,如果很小则可能是小RNA的片段(一部分)。

2.4 Counts file

计数文件,可以用来做差异表达分析。

2.5 MIRNAs directory

包含鉴定出来的miRNA基因簇信息,依据的是“Results.txt”中MIRNA这一列为“Y”(yes)的记录会被确定为一个"new miRNA family"。

每一个文件具体信息包括:这段位置的序列,以点和括号表示的二级结构,miRNA and miRNA-star的位置;下方的记录表示其它的小RNA,当比对到相反链时,用"<"表示;l: length of RNA, a: number of alignments.

2.6 gff3文件

根据Results.txt中DicerCalls这一列是否为'N',得到ShortStack_D.gff3,ShortStack_N.gff3,ShortStack_ALL.gff3这几个粗略注释文件

3. ShortStack详细使用指南

Usage: ShortStack [options] {--readfile <r> | {--bamfile <b> | --cramfile <c>}} --genomefile <g>
<r> : *.fasta or *.fa; colorspace-fasta (.csfasta); *.fastq or *.fq; their gzip-compressed versions(多文件时,空格分隔)
<b> : *.bam
<c> : *.cram
<g> : FASTA formatted (.fa or .fasta) genome file

Alignment Options:
--adapter:至少由8个ATGC碱基组成的序列
--bowtie_cores:bowtie比对时所用线程数
--sort_mem:samtools sort所用的内存
--mismatches:比对时能容忍的错配数,默认是1
--mmap:处理多比对reads的方式,n (none), r (random), u (unique-seeded guide), or f (fractional-seeded guide). default: u
'n,忽略多比对reads;r,随机给多比对reads分配一个位置;
u,假设一条READS比对到了基因组的A, B, C三个位置,在这三个位置上划定50nt的窗口,计算窗口内特异性比对reads的数量。READS更有可能落在数量大的那个窗口对应的基因组位置上;
f,在上面比较,不同的是,除了计算窗口内特异性比对reads的数量,还要计算多比对reads的数量;'
—— 配图1
--bowtie_m:多比对最多比对多少次,默认50;超过此值就定为“未比对上”
--ranmax:当--mmap使用u, f时,也可能出现不能准确定位的情况,比如,同时存在多个权重最大的位置。这时,则会随机分配可能位置,但是也有条件,就是可能的位置数不能大于--ranmax,一旦超过就会定为“未比对上”,默认是3
--align_only:跑完比对+定相(处理多比对reads的方法)之后就停止
--show_secondaries:显示primary alignments以及secondary alignments,如果已经比对定相了,这个选项还有用吗?
Analysis Options:
--dicermin:Dicer酶切割得到的小RNA的最小长度,默认20
--dicermax:Dicer酶切割得到的小RNA的最大长度,默认24 ——配图2
--foldsize:在基因组上搜寻MIRNA基因时,考虑到空间结构的折叠(反向重复,茎环结构),片段大小不应该大于一个阈值,一旦大于则不会被看作是MIRNA,默认是300
--locifile:制表符分隔的普通文本文件。第一列Chr:start-stop,第二列为名称(可选)
--locus:与--locifile作用相似,且不能同时使用。Chr:start-stop[,Chr:start-stop,Chr:start-stop......]
--nohp:不执行MIRNA搜索
--pad:如果小RNA的基因簇之间的距离小于等于pad值则会被合并,默认值75 —— 背景知识见下文引用
--mincov:小RNA的基因簇最少的比对reads数,也可以是标准化值,如rpm, rpmm(这时要带小数点),默认值是0.5rpm
--strand_cutoff:对于基因组上一个区段,其实有两条链,需要根据--strand_cutoff这个阈值来定义正负(+-)链。该区段上超过80%的reads比对到参考基因组链,则称为+链;小于等于20%reads比对到参考基因组链定义为-链
1. 比对算法
2. 来自中国大学MOOC,《基因组学》——复旦大学

关于基因簇:大多数miRNA 基因以单拷贝、多拷贝或基因簇的形式存在于基因组中。基因簇(gene cluster)指基因家族中的各成员紧密成簇排列成大串的重复单位,位于染色体的特殊区域。基因簇少则可以是由重复产生的两个相邻相关基因所组成,多则可以是几百个相同基因串联排列而成。他们属于同一个祖先的基因扩增产物。也有一些基因家族的成员在染色体上排列并不紧密,中间还含有一些无关序列。但总体是分布在染色体上相对集中的区域。——百度词条

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

推荐阅读更多精彩内容