生信技能23 - Samtools提取BWA比对的mapped reads

最新推荐文章于 2025-05-06 17:20:31 发布
最新推荐文章于 2025-05-06 17:20:31 发布
阅读量979 收藏
点赞数
CC 4.0 BY-SA版权
分类专栏： 生信分析项目实战技能合集 文章标签： 数据挖掘 数据分析 bash
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/LittleComputerRobot/article/details/130311284
本文介绍了如何通过Samtools从BWA比对结果中提取mapped reads，首先构建参考序列索引，然后进行bwa比对并转换SAM为BAM格式，最后用samtools提取比对成功的reads。

本文目的：筛选出能比对到参考序列的mapped reads，形成新的fastq文件

参考序列构建索引

bwa index ref.fasta

运行成功后会生成5个结果文件：ref.fasta.amb、ref.fasta.ann、ref.fasta.bwt、ref.fasta.pac、ref.fasta.sa

bwa比对+samtools进行sam2bam转换

利用samtools view 命令可将bwa比对生成的为sam（sequence Alignment mapping)文件，将SAM转换为二进制对应的BAM格式。

bwa mem ref.fasta sample_R1.fq.gz sample_R2.fq.gz | samtools view
了解本专栏 订阅专栏 解锁全文 超级会员免费看
samtools提取bam paired-end reads
qq_21478261的博客
09-15 666
通过指定samtools view 的-f /-F/-G参数的FLAGS即可愉快的从bam中提取预期的reads -f INT only include reads with all of the FLAGs in INT present [0] -F INT only include reads with none of the FLAGS in INT present [0] -G INT only EXCLUDE reads with all of the FLAGs in I
生信技能46 - Call人类线粒体变异和提取chrM变异位点
LittleComputerRobot的博客
05-25 379
使用bwa将样本fastq文件比对到线粒体参考基因组， 并使用bcftools进行call变异。
【bioinfo】根据sam文件中的MD标签判断reads比对情况
lillianqdzpw
10-20 840
MD标签是值：比对上的相对位置信息和错误信息。根据MD标签判断read的比对情况。
anchorwave进行复杂基因组比对（3）
m0_69603624的博客
08-16 1366
anchorwave以玉米B73和Mo17为例对存在倒位的基因组 进行全基因组比对。与其他主流软件相比，具有较好的敏感性和准确率。
BWA比对及Samtools提取目标序列
weixin_56827666的博客
05-13 6523
今天想看一下自己的序列里面会不会有某细菌基因组存在，主要使用BWA和Samtools：bwa主要用于将低差异度的短序列与参考基因组进行比对。主要包含三种比对算法：backtrack、SW和MEM，第一种只支持短序列比对（<100bp），后两种支持长序列比对(70bp~1M)，并支持分割比对（split alignment）。MEM算法是最新的也是官方推荐的。BWA-MEM 是一种新的比对算法，用于将测序 reads 或者组装后 contigs 比对至大型参考基因组，例如人参考基因组。它会自动选择局部比
怎么从bam文件中提取出比对OR没比对上的paired reads | bamToFastq | STAR
weixin_34346099的博客
03-27 2634
折腾这么多都是白瞎，STAR就有输出没有别对上的pair-end reads的功能   参见：How To Filter Mapped Reads With Samtools I had the same issue but with Paired End Reads, and I solved using samtools and bamToFastq. You can find bamT...
samtools提取指定位置核酸或氨基酸序列
Curry_chenhu的博客
08-22 2470
提取序列： samtools faidx input.fa chr1 > chr1.fa samtools faidx input.fa chr1:100-200 > chr1.fa
生信软件46 - 三代测序低深度全基因组测序结构变异SV检测工具NanoVar
LittleComputerRobot的博客
05-06 854
三代测序低深度全基因组测序结构变异SV检测工具NanoVar
生信分析进阶5 - 全外显子组变异检测和ANNOVAR注释Snakemake分析流程
LittleComputerRobot的博客
06-16 740
全外显子组变异检测和ANNOVAR注释Snakemake分析流程
ChIP-seq数据处理流程(附赠长达5小时的视频指导)
生信小博士的博客
12-15 6443
本次给学徒讲解的文章是 : Brookes, E. et al. Polycomb associates genome-wide with a specific RNA polymerase II variant, and regulates metabolic genes in ESCs. Cell Stem Cell 10, 157–170 (2012). 查看文章发现数据是： Polycomb associates genome-wide with a specific RNA polymerase
samtools 序列比对率计算（samtools flagstat）
wt141643的博客
04-29 7647
准备序列比对后生成的 bam 文件或者 .sam 文件 samtools flagstat .bam文件 > flagstat.tax 结果解释 从第一行至第十一行分别表示： QC pass的reads的数量为2499971，未通过QC的reads数量为0，意味着一共有2499971条reads； 重复reads的数量，QC pass和failed 比对到参考基因组上的read...
linux批量截取一个基因上下游序列,利用samtools截取基因组上指定位置的序列
weixin_39793576的博客
05-15 2986
samtools faidx 能够对fasta 序列建立一个后缀为.fai 的文件，根据这个.fai 文件和原始的fastsa文件， 能够快速的提取任意区域的序列用法：samtools faidx input.fa该命令对输入的fasta序列有一定要求：对于每条序列，除了最后一行外， 其他行的长度必须相同,>oneATGCATGCATGCATGCATGCATGCATGCATGCATGCATG...
不会编程，如何快速提取序列
xxxie_的博客
07-25 1548
提取序列是生物信息分析中常见的一个操作，也是学习生物信息编程的入门操作。通常是给定基因ID，然后从一个大的数据集里面提取出匹配ID的序列，包含匹配的序列ID和序列信息，类似于Excel中的Vlookup，但是这里需要一个包含序列ID的列表以及一个包含序列的fasta格式文件。如果不会编程该如何提取呢，今天我们就介绍一些方法。 例如这里有五条序列，我们需要根据基因ID，提取出gene3和gene5的...
bwa序列比对
qq_27390023的博客
10-06 5974
BWA主要用于把基因组测序数据比对到参考基因组上，由三种类似算法组成：BWA-backtrack，BWA-SW和BWA-MEM。首推BWA-MEM。 1. 建立索引 bwa index hg38.fa 注：在同一目录下生成以hg38.fa为前缀的5个文件，输入不能是.gz压缩文件。索引建立好了，bwa比对时不再需要原始的hg38.fa 文件。 2.用默认参数比对 # SE （single end） bwa mem hg38.fa test.fq.gz > test_se.sam bw.
240124学习记录：查看/提取bam文件的reads
dstyaaaa123456的博客
01-24 4254
提取bam文件的reads数量
读取Unique reads
huyongfeijoe的专栏
06-24 2955
从贴好的bam文件中读取unique mapped reads，来进行下一步的分析。 tophat贴的结果对于非unique的reads会反复出现在bam/sam中，因此不挑unique reads后续结果会不准确。 Tophat结果中Unique reads的Quality一栏为255，而BWA贴出的结果此栏会根据reads本身的Phred值变化，因此确定unique的唯一标准是option
宏基因组实战7. bwa序列比对, samtools查看, bedtools丰度统计
热门推荐
刘永鑫的博客——宏基因组公众号
11-13 1万+
前情提要如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章 宏基因组分析理论教程 微生物组入门圣经+宏基因组分析实操课程 1背景知识-Shell入门与本地blast实战 2数据质控fastqc, Trimmomatic, MultiQC, khmer 3组装拼接MEGAHIT和评估quast 4基因注释Prokka 5基于Kmer比较数据集sourmash 5基于Kmer比
samtools常用命令详解
weixin_33883178的博客
03-05 1466
samtools的说明文档：http://samtools.sourceforge.net/samtools.shtml samtools是一个用于操作sam和bam文件的工具合集。包含有许多命令。以下是常用命令的介绍 1. view view命令的主要功能是：将sam文件转换成bam文件；然后对bam文件进行各种操作，比如数据的排序(不属于本命令的功能)和提取(这些操作 是对bam文...
生信基础概念之unique reads VS multi-mapping reads
最新发布
09-02
在生信领域，unique reads（唯一比对的读段）和 multi - mapping reads（多重比对的读段）是在将测序得到的短序列（reads）比对到参考基因组时出现的两种不同情况。 Unique reads 指的是那些在参考基因组中只有一个最佳比对位置的 reads。这些 reads 能够明确地定位到基因组的某一个特定区域，因此在后续的分析中，它们所携带的信息通常比较可靠，能够为基因表达定量、变异检测等分析提供准确的数据支持。例如在基因表达定量分析中，unique reads 可以准确地反映某个基因的转录本丰度。 Multi - mapping reads 则是那些可以比对到参考基因组中多个不同位置的 reads。这种情况通常是由于基因组中存在重复序列，或者基因家族成员之间具有高度的序列相似性导致的。由于这些 reads 无法明确其来源，在后续分析中使用时会带来一定的困难。例如在基因表达定量分析中，很难确定这些 reads 究竟是来自哪个基因。 在代码示例中，以 BWA（Burrows - Wheeler Aligner）比对工具的输出为例，SAM（Sequence Alignment/Map）格式文件中，对于 unique reads 和 multi - mapping reads 有不同的标记。以下是简单的 Python 代码示例，用于统计 unique reads 和 multi - mapping reads 的数量： ```python unique_count = 0 multi_count = 0 with open('alignment.sam', 'r') as f: for line in f: if line.startswith('@'): continue fields = line.strip().split('\t') # 这里简单假设 MAPQ 值为 0 表示 multi - mapping reads if int(fields[4]) == 0: multi_count += 1 else: unique_count += 1 print(f"Unique reads count: {unique_count}") print(f"Multi - mapping reads count: {multi_count}") ```
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

生信与基因组学

每一份鼓励是我坚持下去动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值