Skip to content

Latest commit

 

History

History
187 lines (116 loc) · 12.8 KB

How2Deal.md

File metadata and controls

187 lines (116 loc) · 12.8 KB

解题思路提示

目录

入门题 目录

1. RNA序列翻译成蛋白质 目录

2. 获得反向互补序列 目录

3. 根据序列ID，提取目标序列 目录

（1）解析Fasta/Fastq文件

将给定的Fasta/Fastq文件逐行读入，解析ID行和序列行，若是Fastq文件，还需要解析质量行（它在Fastq文件的第4行）

注意：在Fastq文件中，每条记录的第3行为保留行，均为“+”字符，不包含任何信息，因此可以不需要进行解析，也不需要特意定义变量来保存这一行的信息

解析好之后，以哈希形式（在Perl中对键-值形式数据结构的叫法，在Python中成为字典）来存储解析后的数据，其中key为序列ID，value为序列字符串（若为Fastq文件，需要创建两个哈希，一个哈希保存：ID -> Seq，另一个哈希保存：ID -> Qual）

（2）根据ID list提取目标序列的记录

将给定的ID list文件逐行读入，若该ID所代表的序列在给定的Fasta/Fastq文件中存在，则在前面所创建的哈希变量中就能找到和它对应的键名(key)，则将其写出，否则给出提示：No such sequence!

示例代码，点这里

4. 双端未匹配数据的重新匹配 目录

(1) 解析Fastq文件

对双端的两个Fastq文件分别进行解析，保存成两组哈希，每组哈希有两个哈希变量，分别保存：ID -> Seq 和 ID -> Qual

解析方法与 3. 根据序列ID，提取目标序列中的相同

(2) 找出双端能匹配的fragments将它们分别在输出的两个Fastq文件中在同一行写出

以两组哈希中的任意一组中的一个哈希的键进行遍历，例如以组一的ID -> Seq的键进行遍历，然后在另一组的ID -> Seq中查到对应的键是否存在，若存在，就说明双端是配对的，则进行输出

示例代码，点这里

5. 将输入的大Fasta文件拆分成若干个小Fasta文件 目录

该任务允许用户指定两种拆分方式：指定子文件中的序列条数（为了方便下面的讨论说明，设为N）或者指定最终的输出子文件数（为了方便下面的讨论说明，设为M）——对于用户可能进行的不规范操作，即同时指定这两个条件，很可能这两个条件是相互矛盾的，此时只遵循前一个设定

不论用户最终实际上使用了哪种设定，都需要先解决一个问题：每一个子文件的序列条数n应该是多少条？

  • 如果用户已经直接指定了子文件的序列条数，那么就直接获取，不需要再去求解了，此时n=N；

  • 如果用户指定的是最终生成的子文件数，那么可以通过计算：N=输入Fasta文件的中序列条数 / 子文件数= 输入Fasta文件的中序列条数 / M，因为N可能不是整数，而序列条数是整数，肯定需要对N进行取整，那么如何取整呢？

    输入Fasta文件的中序列条数，可以通过统计输入文件中以'>'起始的行数来确定

    考虑一下上面可采用的商的取整的方案（向上取整ceil or 向下取整floor），以及不同的取整方案可能带来的影响：

    • 若商不为整数，此时若向下取整，则意味着，实际上子文件中的序列条数比期望中的少一点，即n<N，则最终实际拆出的子文件会比用户指定的多至少一个，即n<N，m>=M+1；

    • 若商不为整数，此时若向上取整，则意味着，实际上子文件中的序列条数比期望中的多一点，即n>N，那么最终实际上拆出的子文件数会小于等于用户指定的数量M，即m<=M

      到底最后是等于还是小于，取决于之前的M-1个子文件中总共多写的序列条数是否大于n，即是否满足 (N-n)(M-1) >= n？当(N-n)(M-1) < n时，能保证正好拆出M个子文件，其中前M-1个子文件平均都有n条子序列，最后一个子文件的序列条数不满n条

      通过求解上面的不等式，可以算出n应该满足：n > N(1- 1/M)，而上面的向上取整已经限定了n>N，它满足n > N(1- 1/M)这个要求，所以向上取整能够保证得到的子文件数正好等于M

通过上面的推导论证，每一个子文件的序列条数n应该等于用户指定的N，或者N=输入Fasta文件的中序列条数 / M的向上取整

则剩下要做的就是逐行读入输入Fasta文件，通过对行首起始的>字符的识别来计数当前的Fasta序列的条数i，若i是能被n整除，则说明当前子文件刚好写满，此时要新起一个子文件

示例代码，点这里

进阶题 目录

1. 从Fastq文件中随机抽样一定量的数据 目录

3. 将若干个单样本的表达定量结果汇总成一个大矩阵，即expression profile matrix 目录

本题有两个推荐的实现思路：

  1. 通过构造双重哈希实现

    （1）根据指定的文件夹（下面记作dir）和文件后缀（下面记作pattern），将文件路径为dir/*pattern的文件逐一读入，然后保存为双重哈希形式，即

    %H={feature → {sample → quant}}

    同时记录下样本名的列表 @S，最后输出的矩阵的列所表示的样本的顺序，由样本名列表 @S 确定

    （2）然后，根据上面构造出来的双重哈希 %H 和样本名列表 @S，对每个 feature 逐行写出到输出文件中，若当前 feature 为 i，遍历样本名列表 @S，若当前样本名为 j，则 feature i 在样本 i 的取值记为 nij，缺失值用0填充

    nij = defined(H{i}{j}) ? H{i}{j} : 0

    示例代码：点这里（基于perl的实现）

  2. 先构造初始矩阵，然后再进行填充

    （1）根据指定的文件夹（下面记作dir）和文件后缀（下面记作pattern），将文件路径为dir/*pattern的文件逐一读入，得到 unique feature list 和 sample list，为了方便后面的说明，分别设为变量 @feature_list 和 @sample_list

    （2）根据 @feature_list 和 @sample_list，构造初始矩阵，行数为 @feature_list 的长度，列数为 @sample_list 的长度，矩阵中的每一个元素的值都赋为0（以0为缺省值）

    矩阵的行与 @feature_list 对应，矩阵的列与 @sample_list 对应

    （3）再根据指定的文件夹（下面记作dir）和文件后缀（下面记作pattern），将文件路径为dir/*pattern的文件逐一读入：

    • 根据读入的文件的文件名，获取样本id，将它与 @sample_list 比较，从而确定对应的矩阵的列索引 j；

    • 根据读入文件中第一列的GeneId，将它与 @feature_list 比较，从而确定对应的矩阵的行索引 i；

    获得行索引 i 和列索引 j 之后，就修改矩阵中对应元素的值了

    注意：该方法适用于数据量比较小的情况，当数据量比较大时，推荐用第一种方法实现

    示例代码：点这里（基于R的实现） 点这里（基于python的实现，该示例脚本由Zhongyi Hua同学提供）

4. 利用Needleman–Wunsch 算法来编写一个简单的全局比对程序目录

（1）熟悉Needleman–Wunsch算法，可以在纸上利用二维矩阵画出最优的比对路线。

（2）然后你在纸上如何画的，你就如何写程序实现。如，先初始化一个二维矩阵，填写必要数据，按照公式，将矩阵中的每一小格都填好数据，其中每一小格需记录好，累积分数和最优来源。数据填完就开始回溯找出比对方案.

示例代码：点这里

挑战题 目录

4. 相似数组搜索 目录

算法实现的逻辑为：

对应一个给定的查询数组（来自于S.txt），逐一对每个数值进行比较；

对应当前目标数值x，快速在P.txt文件中，快速找到落在 [x-d, x+d ] (其中d表示定义的最远相似距离) 的数值集合 N={n1, n2, …}，以及集合N中每个数值对应的P.txt中数值来源A={a1, a2, …}， 然后对找到的数值计分，计一分；

当一个数组遍历比较结束后，找出得分最高的P.txt中的数组

那么这个问题的关键就在于：怎么快速找到落在[x-d, x+d ] 的数值集合N={n1, n2, …}，以及集合N中每个数值对应的P.txt中数值来源A={a1, a2, …} ？

关键在于对 P.txt 的高效索引，主要采用了倒排索引 + 哈希链表的方法

  1. 采用了数字到数组的倒排索引

    对 P.txt 文件中的每个数值，建立数值指向其来源数组的倒排索引，一旦定位到了这个数值就可以马上得到这个数值来源于哪个数组

  2. 以哈希链表形式来组织P.txt中的数值（即倒排索引的字典）

    目标搜索的数值范围是 [x-d, x+d]，但是一开始的搜索起点 start 肯定要 start <=x-d ，最后的搜索终点end肯定要保证 end >= x+d，即实际的搜索范围 [start, end] 肯定要包含 [x-d, x+d]，而为了提高搜索效率，要尽量缩小搜索范围，减少不必要的搜索，对此我采用哈希链表的方法：

    以d的从左到右第一位非0有效数字的更高一位作为一个粗精度，例如d=0.01，则我们以0.1作为一个粗精度，将落在[X.1, X.2 )
范围的数值存在一个递增链表里，并用X.1这个值作为哈希的键，指向这个递增链表，则根据键，就可以知道其所存储的链表的数值区间

对于一个给定数值x，它的目标搜索区间是[x-d, x+d]，这个区间要么落在一个咱们上面定义的区间内，要么横跨相邻的两个区间，
则在一个或两个区间内搜索就可以了

以上是算法性能提升的主要实现原理，除此之外，还在以下几个小细节处，进行了优化：

  • 得到最高得分的数组

    当对一个数组中的每个数值进行相似比较计分后，最后需要从计分结果里把得分最高的那个数组找出来，一种最简单粗暴的方法是：最后遍历整个得分结果，找出最高分的那个数组

    但是这样做太低效，我在这里设置了一个变量来记录进行每一个数值的比较打分后，当前得分最高的数组，若有更高分，则对它进行更新

    这样一个查询数组里的每个数值比较打分结束后，直接就得到了得分最高的数组

示例代码：点这里

补充学习资料：