处理rilseq实验结果

RILseq的Python项目详细描述


=====
rilseq
======
intention
----
这个包可以用来分析rilseq实验。它是为一个原核基因组而写的,没有剪接连接图和一些附加的特征。rilseq在melamed等人,分子细胞63(2016),第884-897页(http://www.cell.com/molecular cell/fulltext/s1097-2765(16)30413-0)中有描述。


该软件包处理将fastq文件处理到相互作用的成对rna和一些统计数据的不同阶段。它*不*处理质量问题、适配器移除等。因此,在应用此软件包之前,应使用cutadapt或等效软件处理fastq文件。

在这个自述文件中,假设读取是成对结束的,但是所有的命令都可以用于单端映射。


run with::

$map_single_fragments.py genome.fa-g annotations.gff-1 lib1_1.fastq[.gz]lib2_1.fastq[.gz]lib3_1.fastq[.gz]-2 lib1_2.fastq[.gz]lib2_2.fastq[.gz]-d output_dir-o output_head[-r]-m max_mismatches

此脚本使用bwa将读取映射到基因组。可以为每个库指定两个文件,也可以使用-2或仅使用-1指定一个文件。-2中库的顺序应与-1中的顺序相同,但如果某些库是单端排序的,则顺序可能更短(在本例中,单端应该是-1列表中的最后一个)。还有其他参数可以更改,有关其他帮助,请参见-h。



chimera mapping
----
在为每个库创建BAM文件(使用map_single_fragments.py或任何其他mapper生成)之后,我们可以查找嵌合体读取。

文件。在成对末端测序的情况下,取前25个核苷酸(或按-l参数指定)。在单端测序的情况下,每次读取的前25个和最后25个,确保两个区域不重叠,否则不会得到任何结果。

bam文件中的所有读取都映射到基因组并写入结果文件,除非定义了-a文件名或-a,在第一种情况下,单个片段将写入指定的文件。单个片段用于统计测试,并将标记为"单个",因此将根据它们测试交互作用,但不会使用Fisher精确测试。由于某种原因,如果一个片段是一个单一的一致映射,而这里的单一映射不是一个选项,那么它将被忽略。


25 nt长的序列使用滤尘器进行筛选,以低复杂度去除读操作。灰尘过滤器的默认阈值为10,当0时,可以使用--dust_thr参数更改该阈值。提取两端后,使用BWA将它们映射到基因组,并再次筛选,以查看它们是否可以在同一转录本上。为了做到这一点,我们允许相对大量的不匹配(默认情况下为3,设置为--max_mismatches),并测试每个读取映射到的位置的任何组合是否可以从同一个转录本中产生。我们移除相距1000 nt的成对读取,并按预期或相反顺序映射到相反的方向(循环rnas产生的读取,使用--keep_circular省略此选项)。如果给定-t参数,则测试这些对,看它们是否可能来自同一个转录本,即使它们的距离大于1000nt。此选项在筛选有时来自长转录本的rrna时非常有用。

如果读取的映射至多有1个不匹配项(设置为允许的不匹配项),则收集每个读取的染色体上的最低位置。

e和单词chimera/single根据片段的性质而定。作为输入提供的所有BAM文件都将连接到同一个输出文件,它们可以使用读取名称进一步分离。或者,您可以单独运行每个BAM文件,并对所有读取的文件进行分类。

运行时使用:


$map_chimeric_fragments.py genome.fa lib1_bwa.bam lib2_bwa.bam…

由于测序结果包含非特异性嵌合体,因此需要另一个阶段去除实验噪声。这是通过比较支持交互的读取数与随机期望的读取数来实现的。简单地说,我们计算一个2x2列联表,其读取次数如下:



======================================================================================================================================================n
==
==


l和m还包括单个读取的数量,统计测试两个区域的交互频率是否比随机预期的更频繁。
如果k大于预期,则这两个区域可能在实际交互中
*在活体内*。比值比用(k/l)/(m/n)计算,p值用fisher精确检验计算,仅当k大于预期时进行检测(单尾检测)

它们(设置为-min_int)。如果p值小于0.05(-max_pv),则报告此对。为了避免重新报告相邻区域,这些区域可以扩展到500个NTS(-maxsegs),并且将报告具有较低p值的区域的组合(在Bonferroni校正之后)。

相对于细胞内所有rna而言,hfq较小,这意味着相互作用将没有或几乎没有影响。为了纠正这个问题,如果给出总rna实验的结果,就可以计算出hfq限制的分数的近似值。使用--总RNA后跟以逗号分隔的索引BAM文件列表。脚本将计算每个区域总rna的读取次数,并估计hfq的丰度,比值比将乘以这两个RNAs的值,并在"总标准化比值比"栏中报告。


在确定区域后,用区域内的基因等其他数据修饰它们,如果这是一个已知的靶点,则这两个基因等。

-为了得到基因注释,你应该得到你生物体的ecocyc平面文件,它们需要注册,用-ec-dir指向数据目录。染色体的名称可能与bam文件(用于映射的genome.fa文件)和ecocyc文件不同。您可以使用逗号分隔的名称列表(其中ecocyc中的名称跟在bam文件中的名称后面)将一个字典从bam提供给ecocyc脚本。

除了打印交互,此脚本还可以使用rnaup计算交互自由能(仅限版本1,版本2不起作用)如果--shuffles是>;0,则使用shuffled序列计算此能量的p值。



生成绘图和轨迹
——
脚本plot_circos_plot读取map_chimeric_fragments.py的输出到
生成染色体上各区域之间的相互作用列表。它不能
显示两条染色体之间的相互作用。

连同data/e戋coli戋u k12目录中的conf文件和该目录中的短脚本plot戋u interactions.sh,您可以绘制与srna之间的相互作用,基因组上的rrnas和trnas。

运行:

plot_interactions.sh interactions.txt interactions_plot.png

(其他格式也可用,如SVG)

可以使用generate_bed_file_of_endpoints.py将嵌合体片段的读取写入bed文件。该文件打印BAM文件中每个读取的位置,该文件被发现是嵌合的。有一个选项可以只打印作为重要交互的一部分的片段,使用-s interactions_file.txt来完成。当使用-s时,您可以指定一个基因名(ecocyc id),并生成一个bed文件,其中一个片段被映射到该基因(-e id)。运行generate_bed_file_of_endpoints.py-h以获取完整的文档。


数据文件
——
此软件包适用于大肠杆菌k12(refseq nc_000913.2 genome和refseq nc_000913.3 genome)。data
目录包含两个单独的子目录,分别称为ver2和ver3,用于两个基因组版本,其中
包括genome\*.fa、ecocyc genes gff文件和ecocyc transcripts gff文件。这些文件以及ver2和ver3目录中的其他文件分别基于ecocyc版本19.0和20.0,其中包括根据sri international许可证从biocyc(tm)途径/基因组数据库获得的数据。
图片:biocyc logo.jpg
在使用之前,应该使用bwa index genome.fa索引基因组。这两个gff文件可以使用脚本生成::

generate_transcripts_gff.py ecocyc_data_dir

和::


genes_gff.py ecocyc_data_dir

ver2 data目录中还有两个额外的文件:从ecocyc获取的目标的精确列表
稍有改动和一个rep元素表(用于结果注释),此表下载自:http://ecocyc.org/group?id=biocyc14-8223-3640227683

scipy
-biopython

该项目托管在github上:https://github.com/asafpr/rilseq

欢迎加入QQ群-->: 979659372 Python中文网_新手群

推荐PyPI第三方库


热门话题
java SUN次要代码309含义   java避免为空元素生成XML自关闭标记,并生成自定义的<XML>开始标记   java使用json和restful将数组数据从本地sqlite数据库插入SQL Server   java Spring Boot 1.5.9字符编码问题   LInkedIn讨论中的java 401错误   位图Java:检查多个位向量/位集是否相交的最快方法?   macos如何让Java应用程序以图标出现在Mac OS X dock中   java如何删除netbeans中的@SuppressWarnings(“未使用的”)?   apachestorm中的java自定义序列化   java可以退出代码还是应该终止main?   递归如何在Java中递归地绘制简单的线条?   unicode在Java中确定特定字体是否可以呈现特定字符   打开并阅读带有Selenium/Katalon(Java)特定标题的电子邮件文本(来自Gmail)