使用hisat2开始分析

Posted 2023-04-30

参考技术A 1.构建索引

$HISAT2_HOME/hisat2-build $HISAT2_HOME/example/reference/22_20-21M.fa \

--snp $HISAT2_HOME/example/reference/22_20-21M.snp 22_20-21M_snp

2.序列比对

a 单端序列

$HISAT2_HOME/hisat2 -f \

-x $HISAT2_HOME/example/index/22_20-21M_snp \

-U HISAT2_HOME/example/reads/reads_1.fa \

-S eg1.sam

b 双端序列

$HISAT2_HOME/hisat2 -f \

-x $HISAT2_HOME/example/index/22_20-21M_snp \

-1 $HISAT2_HOME/example/reads/reads_1.fa \

-2 $HISAT2_HOME/example/reads/reads_2.fa \

-S eg2.sam

3. 使用 SAMtools/BCFtools 进行下游分析

SAMtools是一组用于操作和分析SAM和BAM比对文件的工具。BCFtools是一组用于calls变异和操作VCF和BCF文件的工具，它通常与SAMtools一起发布。同时使用这些工具允许您操作从SAM格式的比对到VCF格式的calls变异。本示例假设已经安装了samtools和bcftools，并且包含这些二进制文件的目录位于PATH环境变量中。

a  双端序列比对

$HISAT2_HOME/hisat -f \

-x $HISAT2_HOME/example/index/22_20-21M_snp \

-1 $HISAT2_HOME/example/reads/reads_1.fa \

-2 $HISAT2_HOME/example/reads/reads_2.fa \

-S eg2.sam

b SAM文件转化为BAM文件

samtools view -bS eg2.sam > eg2.bam

c  将原始 BAM 文件转化为排序过的 sorted BAM 文件

samtools sort eg2.bam -o eg2.sorted.bam

d  我们现在有一个已排序的 BAM 文件，名为 eg2.sord.bam 。bam文件是一种有用的格式，因为比对内容是 (a) 压缩的，这便于长期存储， (b)是 排过序的，这便于 calls 变异 。要生成 calls 变异的 VCF 格式文件，请运行 :

samtools mpileup -uf $HISAT2_HOME/example/reference/22_20-21M.fa eg2.sorted.bam | bcftools view -bvcg - > eg2.raw.bcf

e  要查看 bcf 格式文件

bcftools view eg2.raw.bcf

从动态规划到配对序列联配

如何找到一个物种的同源基因呢？通常而言，我们都会以BLAST起手，根据序列相似度找到目标物种可能的部位。

基因组重测序得到二代短读长序列(read)如何鉴定变异变异？通常而言，我们会用BWA/Bowtie2起手，将100到150bp长度的read回帖到参考基因组，然后分析测序结果和参考基因组有哪些不同。

转录组测序得到二代短读长序列如何定量呢？通常而言，我们会以STAR/HISAT2起手，将100到150bp长度的read回帖到基因组上，根据基因注释信息统计每个基因的read数。

PacBio/Nanopore的三代测序如何组装？我们也是通过序列相互比对，把相似的序列堆叠在一起进行纠错，然后再相互比对，根据overlap构建graph，然后对graph进行解析得到组装结果。

综上，我们的研究或多或少都和序列比对有关，依赖于序列比对的准确性。如果我们发现一个软件的比对不是特别好，比如说没有BLAST到我们需要区域，通常我们也不会想自己开发一个软件，而是换一个软件。当然我大部分时候也不会想着自己去造轮子，除非没有合适的轮子（哪个新造轮子的不这样子想呢？）。不过，有些时候你需要真的造轮子才能加深你对某个知识点的理解，而这里我就想用自己学的动态规划知识来写一个配对序列联配小程序。

动态规划

第一次看到动态规划(dynamic programming)的时候，我对这个算法就有了不切实际的幻想。"动态"，"规划"，这两个都是好词啊，从里面我看到了人生的哲学，这是让我们不要以静态的观点看问题，要结合实际情况随机应变，规划好未来呀！

接着，我看了具体的导学问题，求解是斐波那契数列第n项，好简单呀。再看练习题，计算从一个字符串到另一个字符串的最小操作次数（编辑距离），我蒙住了，完美没有思路。这种状态俗称，“一看就会，一做就废”。最重要的是，我发现无论是简单的题目还是难题，其实都没有让我想到动态规划中，动态和规划两个词。

为了解决这个疑惑，我们需要对动态规划有正确的认识。不要对名字过多的遐想，而是从它的定义出发。

Dynamic programming is both a mathematical optimization method and a computer programming method. The method was developed by Richard Bellman in the 1950s and has found applications in numerous fields, from aerospace engineering to economics. In both contexts it refers to simplifying a complicated problem by breaking it down into simpler sub-problems in a recursive manner. While some decision problems cannot be taken apart this way, decisions that span several points in time do often break apart recursively. Likewise, in computer science, if a problem can be solved optimally by breaking it into sub-problems and then recursively finding the optimal solutions to the sub-problems, then it is said to have optimal substructure.

根据维基百科的定义，我们知道动态规划的本质是将复杂问题拆分成简单的子问题，从而简化问题(simplifying a complicated problem by breaking it down into simpler sub-problems )。动态规划并非万能，它只能解决一些存在最优子结构的问题(optimal substructure).

因此，动态规划并没有它的名字那么玄奥，也不是我们想的那么万能。后续为例避免大家看见名字就想入非非，我们一律用简写DP来表示该算法。

DP的关键是求解子问题的时候，能够重复利用(reuse)已经求解的子问题结果，而不是从头计算，因此降低了计算的时间复杂度（但是提高了空间复杂度）。它有两种形式

• 自顶向下（递归+记忆化）

• 自底向上（递推+状态表）

接下来，我们以最经典的斐波那契数列求解来介绍这两种DP。

斐波那契数列

我们都知道斐波那契数列的递推公式为f(n) = f(n-1) + f(n-2), 其中F(0) = 0, f(1) =1. 于是我们就可以通过递归的方式实现，下面是一段R语言的代码

 
   
   
 

  
    
    
  
fib <- function(x){
  
    
    
  
 if (x <= 1) return(x)
  
    
    
  


  
    
    
  
 return( fib(x-1) + fib(x-2))
  
    
    
  
}
 
   
   
 

代码的缺陷在于，越到后面计算量越大，所以计算时间呈现出一种指数增长的状态。

 
   
   
 

  
    
    
  
run_time <- vector(length = 32)
  
    
    
  


  
    
    
  
for (i in seq(32)){
  
    
    
  
 run_time[i] <- system.time(fib(i))[3]
  
    
    
  
}
  
    
    
  


  
    
    
  
plot(x=1:32, run_time)
 
   
   
 


为了找到原因，我们需要对计算过程进行拆解。从递归树展开中，我们可以发现很多计算出现了重复，比如说fib(6)中计算了fib(5), 那么在计算fib(7)的时候，就不需要再次计算fib(5).

为了避免这种情况，我们就可以建立一个"备忘录"来复用已经计算的结果.

 
   
   
 

  
    
    
  
fib2 <- function(x){
  
    
    
  
 if (x <= 1) return(x)
  
    
    
  
 if ( ! is.null(memo[x][[1]]) ){
  
    
    
  
 return(memo[[x]])
  
    
    
  
 }
  
    
    
  
 memo[[x]] <<- fib2(x-1) + fib2(x-2)
  
    
    
  
 return(memo[[x]] )
  
    
    
  
}
 
   
   
 

然后运行它， 你会发现计算时间是一个常数

 
   
   
 

  
    
    
  
run_time <- vector(length = 32)
  
    
    
  


  
    
    
  
memo <- list()
  
    
    
  
for (i in seq(32)){
  
    
    
  
 run_time[i] <- system.time(fib2(i))[3]
  
    
    
  
}
  
    
    
  


  
    
    
  
plot(x=1:32, run_time)
 
   
   
 

这就是自顶向下的DP。

当然，我们也可以选择自底向上的递推，依次向上算出f(0), f(1)..f(n).

 
   
   
 

  
    
    
  
fib3 <- function(x){
  
    
    
  
 if ( x <= 2) return(x);
  
    
    
  
 vec <- vector(length = x)
  
    
    
  
 vec[1] <- 1
  
    
    
  
 vec[2] <- 2
  
    
    
  
 for (i in seq(3,x)){
  
    
    
  
 vec[i] <- vec[i-1] + vec[i-2]
  
    
    
  
 }
  
    
    
  
 return(vec[x-1])
  
    
    
  
}
 
   
   
 

每次计算时间也是常数级别。就不再作图展示。

通常而言，我们会更偏向于自底向上的递推方法，可以避免系统调用栈的额外开销，后续说的DP都会以自底向上的递推方法进行介绍。

小结

这篇文章主要介绍了动态规划定义，以及以一个经典的斐波那契数列求解介绍了两种常见DP形式，下一部分则是用另外一个例子，编辑距离，来介绍DP编程的标准步骤。