使用bwa进行序列比对
001、
bwa mem -t 4 -k 32 -M -R "@RG\tID:name\tSM:name\tPL:illumina\tLB:name\tPU:name" reference.fna sm.clean.1.fastq.gz sm.clean.2.fastq.gz | samtools view -Sb - > sm.bam
mem:mem比对算法
-t:指定线程数
-k: (这个参数可以不设置)最小的种子长度(minimum seed length),bwa采用seed-and-extend的策略,在seed阶段,bwa取read的碱基片段在reference上进行精确匹配,并选择满足一定长度要求和匹配次数的read片段作为seed,这个阶段算法的核心是基于FM-index的精确匹配; 在extend阶段,bwa利用smith-waterman算法将seed在read和reference上向两边延伸比对(容忍gap),进而找到整个read在reference上符合条件的全局匹配。
-M: bwa mem在比对的时候, 对于一条序列同时比对到基因组不同区域的情况,bwa认为都是最优匹配,但是会和picard不兼容,影响后面gatk的检测,这个时候可以设置-M选项,将较短的比对标记为次优,与picard兼容。
-R:Read Group的字符串信息,这是一个非常重要的信息。
(1) ID,这是Read Group的分组ID,一般设置为测序的lane ID(不同lane之间的测序过程认为是独立的),下机数据中我们都能看到这个信息的,一般都是包含在fastq的文件名中;
(2) PL,指的是所用的测序平台,这个信息不要随便写!特别是当我们需要使用GATK进行后续分析的时候,更是如此!这是一个很多新手都容易忽视的一个地方,在GATK中,PL只允许被设置为:ILLUMINA,SLX,SOLEXA,SOLID,454,LS454,COMPLETE,PACBIO,IONTORRENT,CAPILLARY,HELICOS或UNKNOWN这几个信息。基本上就是目前市场上存在着的测序平台,当然,如果实在不知道,那么必须设置为UNKNOWN,名字方面不区分大小写。
(3) SM,样本ID,同样非常重要,有时候我们测序的数据比较多的时候,那么可能会分成多个不同的lane分布测出来,这个时候SM名字就是可以用于区分这些样本。
(4) LB,测序文库的名字,这个重要性稍微低一些,主要也是为了协助区分不同的group而存在。文库名字一般可以在下机的fq文件名中找到,如果上面的lane ID足够用于区分的话,也可以不用设置LB;
除了以上这四个之外,还可以自定义添加其他的信息,不过如无特殊的需要,对于序列比对而言,这4个就足够了。这些信息设置好之后,在RG字符串中要用制表符(\t)将它们分开。
管道后接samtools:
-S:表示读入文件为sam格式
-b:表示输出文件为bam格式
。
参考:https://developer.aliyun.com/article/1242282
。
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 震惊!C++程序真的从main开始吗?99%的程序员都答错了
· 【硬核科普】Trae如何「偷看」你的代码?零基础破解AI编程运行原理
· 单元测试从入门到精通
· 上周热点回顾(3.3-3.9)
· winform 绘制太阳,地球,月球 运作规律
2022-09-22 linux 中统计相同序列出现的次数
2022-09-22 sbatch命令在集群递交任务模板
2022-09-22 linux 中 利用命令向文件的末尾添加空行
2022-09-22 linux中basename命令
2022-09-22 linux 中 date +%s 获取1970年以来的秒数