宏基因组实战7. bwa序列比对, samtools查看, bedtools丰度统计

最新推荐文章于 2024-03-14 11:27:22 发布
最新推荐文章于 2024-03-14 11:27:22 发布
阅读量1.4w 收藏 49
点赞数 6
分类专栏： 宏基因组 宏基因组分析
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/woodcorpse/article/details/78523035
版权

前情提要

如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章

测试数据

刘博士帮助把测试数据建立了一个百度云同步共享文件夹，有非常多的好处，请读完下文再决定是否下载：
1. 下载被墙的数据；很多数据存在google, amazon的部分服务器国内无法直接下载，而服务器一般科学上网不方便，下载数据困难。大家下载失败的数据请到共享目录中查找；
2. 预下载好的软件、数据库；有很多需要下载安装、注册的软件(在线安装包除外)，其实已经在共享目录了，节约小伙伴申请、下载的时间；
3. 数据同步更新；任何笔记或教程不可避免的有些错误、或不完善的地方，后期通过大家的测试反馈问题，我可以对教程进行改进。共享目录不建议全部下载或转存，因为文件体积非常大(目前>18G)，而且还会更新。你转存的只是当前版本的一个备份，就不会再更新了。建议直接在链接中每次逐个下载需要的文件，也对文件有一个认识过程。
4. 方便结果预览和跳过问题步骤；服务器Linux在不同平台和版本下，软件安装和兼容性问题还是很多的，而且用户的权限和经验也会导致某些步骤相关软件无法成功安装(有问题建议选google、再找管理员帮助；想在群里提问或联系作者务必阅读《如何优雅的提问》)。在百度云共享目录中，有每一步的运行结果，读者可以下载查看分析结果，并可基于此结果进一步分析。不要纠结于某一步无法通过，重点是了解整个流程的分析思路。

最后送上本教程使用到的所有文件同步共享文件夹链接：http://pan.baidu.com/s/1hsIjosk 密码：y0tb 。

bwa序列比对

安装bwa

# 工作目录，根据个人情况修改
wd=~/test/metagenome17
cd $wd
sudo apt-get install bwa samtools

下载数据

mkdir mapping
cd mapping
# 可以接之前的学习，或在百度云中下载
cp ../data/*.pe.fq.gz ./ # 引处如果ln硬链解压不允许
# 逐个解压，直接gunzip *.gz批量也不成功
for file in *fq.gz
do
gunzip $file
done

# 无法下载找百度云
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/subset_assembly.fa.gz
gunzip subset_assembly.fa

序列比对

# 建索引
bwa index subset_assembly.fa

# 双端合并序列，使用-p参数比对
for i in *fq
do
# unrecognized command 'mem' 版本过旧，qiime索引了旧版本，bwa改为绝对路
/usr/bin/bwa mem -p subset_assembly.fa $i > ${i}.aln.sam 
done

samtools操作比对结果

samtools可视化比对结果

# 参考序列建索引
samtools faidx subset_assembly.fa

# 压缩sam为bam用于可视化
for i in *.sam
do
 samtools import subset_assembly.fa $i $i.bam
 samtools sort $i.bam -o $i.bam.sorted.bam
 samtools index $i.bam.sorted.bam
done

# 可视化
# 按contig的reads数量排序，找高丰度的查看
grep -v ^@ SRR1976948.abundtrim.subset.pe.fq.aln.sam | cut -f 3 | sort | uniq -c | sort -n | tail
# Pick one e.g. k99_13588.

# 查看k99_13588序列400bp开始
samtools tview SRR1976948.abundtrim.subset.pe.fq.aln.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:400
# 方向可以上下左右移动查看，q退出

# 另一个窗口打开另一个样品，便于比较
samtools tview SRR1977249.abundtrim.subset.pe.fq.aln.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:400

image
在两个终端(Xshell)中用samtools查看不同样品的比对结果序列分析情况

bedtools丰度估计

我们将会用比对结果估计组装序列的覆盖度

使用bedtools比对，常用 http://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html

# 安装程序
sudo apt-get install bedtools

# 使用genomeCoverageBed定量基因
for i in *sorted.bam
    do
    genomeCoverageBed -ibam $i > ${i/.pe*/}.histogram.tab
done

我们看一下结果：

k99_5   6   87  258 0.337209
k99_5   7   18  258 0.0697674
k99_5   8   11  258 0.0426357
  1. Contig name
  2. Depth of coverage 覆盖深度
  3. Number of bases on contig depth equal to column 2
  4. Size of contig (or entire genome) in base pairs
  5. Fraction of bases on contig (or entire genome) with depth equal to column 2

To get an esimate of mean coverage for a contig we sum (Depth of coverage) * (Number of bases on contig) / (Length of the contig). We have a quick script that will do this calculation.

计算平均深度

wget https://raw.githubusercontent.com/ngs-docs/2017-cicese-metagenomics/master/files/calculate-contig-coverage.py
# 安装过可跳过
sudo pip install pandas
for hist in *histogram.tab
do
    python calculate-contig-coverage.py $hist
done

产生结果文件 SRR1976948.abundtrim.subset.histogram.tab.coverage.tab：

k99_100 18.200147167
k99_1000    52.6492890995
k99_10000   4.97881916865
k99_10008   16.7718223583
k99_1001    62.2604006163

可选分析：末质控数据结果如何？

# 下载原始数据
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_1.fastq.gz
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_2.fastq.gz

# 如果有，可复制过来
cp ../data/SRR1976948_?.fastq.gz .

# 解压变只提取前20万行
gunzip -c SRR1976948_1.fastq.gz | head -800000 > SRR1976948.1
gunzip -c SRR1976948_2.fastq.gz | head -800000 > SRR1976948.2

# 比对
bwa aln subset_assembly.fa SRR1976948.1 > SRR1976948_1.untrimmed.sai
bwa aln subset_assembly.fa SRR1976948.2 > SRR1976948_2.untrimmed.sai
bwa sampe subset_assembly.fa SRR1976948_1.untrimmed.sai SRR1976948_2.untrimmed.sai SRR1976948.1 SRR1976948.2 > SRR1976948.untrimmed.sam

# 压缩、排序、索引
i=SRR1976948.untrimmed.sam
samtools import subset_assembly.fa $i $i.bam
samtools sort $i.bam -o $i.bam.sorted.bam
samtools index $i.bam.sorted.bam

# 查看
samtools tview SRR1976948.untrimmed.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:500

比较这些没有trim质控的数据，看看和原来的结果有什么不同？

猜你喜欢

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外七十多位PI，七百多名一线科研人员加入。参与讨论，获得专业指导、问题解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职务”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。
image

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”
image

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA

刘永鑫Adam
关注
  • 6
    点赞
  • 踩
  • 49
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
基于GPU的BWA序列比对算法分析与加速.pdf
09-25
基于GPU的BWA序列比对算法分析与加速.pdf
bwa基因比对程序
03-16
bwa算法，经过bwa算法执行后会得到sam文件，sam文件用于之后的gatk工作
BWA: 快速准确的比对短序列到长序列的工具
最新发布
gitblog_00012的博客
03-14 399
BWA: 快速准确的比对短序列到长序列的工具 BWA 是一个用于将短序列（如 Illumina 测序数据）比对到长序列（如基因组）上的工具。 什么是 BWA？ BWA 是一种基于 Burrows-Wheeler 转换和哈希表的数据结构的比对算法。它可以快速、准确地将短序列比对到长序列上，并且支持多种不同的比对模式，包括全局比对、局部比对和种子比对等。 BWA 可以用来做什么？ BWA 主要用于生物...
测序深度的计算，你真的掌握了吗
庐州月光的博客
03-02 8642
欢迎关注”生信修炼手册”!计算测序深度是数据分析中的一个常规操作，常用的工具有以下几种1. samtools depth命令如下samtools depth input.bam >...
高通量测序分析工具Bedtools使用介绍
悟道西方
04-03 1万+
http://blog.genesino.com/2018/04/bedtools/ Bedtools是处理基因组信息分析的强大工具集合，其主要功能如下： bedtools: flexible tools for genome arithmetic and DNA sequence analysis. usage: bedtools <subcommand> [option...
基于BWA，Bowtie2，Salmon和SAMtools、checkm等工具计算宏基因组学序列分析中Contigs与Genes在样品中的丰度，多种计算方式和脚本对比（20231217更新）
zrc_xiaoguo的博客
12-16 2411
计算contigs的丰度一般使用assembly的结果，计算基因风度时一般用prodigal的结果，prodigal结果中建议同时输出蛋白序列和核酸序列文件，基因注释是一般使用diamond需要使用蛋白序列，而这里计算丰度需要与原始核酸序列进行比对，所以这里要用核酸序列，prodigal输出的核酸序列和蛋白序列id一样，所以只需要最后以序列id进行mapping就可以了。
基因数据处理10之BWA处理速度太慢
Keep Learning
03-13 2738
基因数据处理9之BWA处理速度太慢： hadoop@Mcnode1:~/cloud/adam/xubo/data/data_HDFS/GRCH38/GCA_000001405.15_GRCh38/test20160310$ bwa aln GCA_000001405.15_GRCh38/GCA_000001405.15_GRCh38_full_analysis_set.fna SRR003161
如何获得hg38外显子的bed文件？
weixin_40594350的博客
04-18 3794
一、bed文件介绍 bed文件是一种记录基因组不同（功能）区域在基因组上的位置以及其它注释信息的文本文件。它包含了由空格或者tab分隔的不同列，以记录不同的信息，每一行对应一个区域。它最早出现于人类基因组计划中，后被广泛应用。因为它不直接在基因组上进行标记和修改，在使用上更具效率。 bed文件最开始并没有一个标准的格式，因此 UCSC Genome Browser 对它的描述逐渐成为了大家的参考标准。它最少为3列，最多可为12列。bed文件辅助 UCSC Genome Browser 对不同片段进行可视
生信技能-高通量测序工具bam、samtools、bedtools及conda的下载和安装
weixin_43664814的博客
03-15 2117
用于建立 index；基于 BWT 算法，将 reads 比对到参考基因组；最新版本 bwa-mem2，Intel实验室对计算效率进行了优化。
BWA序列比对方法丨针对较大基因组的并行计算和性能优化方式，利用多线程和负载均衡策略提高效率
青笋的博客
07-02 1756
设置reads标头，也就是sam文件中的RG部分，为什么要设置RG表头呢，因为同一样品可能包括多个测序结果，来自不同lane，不同文库，或者不同样品的比对结果合并到同一个文件中进行处理，就需要通过RG进行标记区分。BWA-MEM默认使用较大的内存来处理测序数据，但在处理大数据量时，可以通过调整BWA的-m参数来限制内存使用。特殊情况下，可以使用 –M 选项来将。BWA-MEM算法是BWA中最常用和最适用于长读段的比对算法，它能够对较长的测序读段执行准确的比对，并具有较好的插入片段和结构变异识别能力。
基因二代测序-分别对blast和bwa比对结果进行统计比较-20230506
12-15
基因二代测序-分别对blast和bwa比对结果进行统计比较-20230506
SqueezeMeta:用于宏基因组学分析的完整管道
04-09
SqueezeMeta：从读取到垃圾箱的全自动宏基因组学流水线 在以下找到SqueezeMeta论文： ://www.frontiersin.org/articles/10.3389/fmicb.2018.03349/full 在以下位置找到有关如何分析SqueezeMeta输出的第二篇论文： :...
BWA 生物信息学工具箱支持包：BWA 将 DNA 序列与参考基因组对齐。-matlab开发
05-31
BWA 是用于将 DNA 序列与参考基因组比对的软件包。 序列比对允许生物学家确定 DNA 片段的基因组来源。 这是分析许多生物信息学数据集的必要第一步，包括 RNA-seq、ChIP-seq 和 DNase-seq。
生信技能23 - Samtools提取BWA比对的mapped reads
LittleComputerRobot的博客
04-22 543
本文目的：筛选出能比对到参考序列的mapped reads，形成新的fastq文件。
宏基因组实战8. 分箱宏基因组binning, MqaxBin, MetaBin, VizBin
刘永鑫的博客——宏基因组公众号
11-21 1万+
前情提要如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章 宏基因组分析理论教程 微生物组入门圣经+宏基因组分析实操课程 1背景知识-Shell入门与本地blast实战 2数据质控fastqc, Trimmomatic, MultiQC, khmer 3组装拼接MEGAHIT和评估quast 4基因注释Prokka 5基于Kmer比较数据集sourmash 5基于Kmer比
Trimmomatic、bowtie2、samtools和bedtools安装过程全记录（已全部安装成功）
qq_42940869的博客
11-08 5902
首先你需要有一个linux的环境 其次在linux系统下用浏览器打开trimmomatic的下载链接，下载最新版本对应的binary下载链接http://www.usadellab.org/cms/index.php?page=trimmomatic下载完成后，将其解压到一个新的文件夹中，然后复制trimmomatic-0.39.jar的路径，打开终端，输入jar <路径>/trimmomatic-0.39.jar， 例如：jar /home/sai/workspace/go/src/pro
BWA MEM比对人类参考基因组详解
睿智如水
11-10 3万+
在获得下机数据后，做的第一步是质控。质控工具有很多，这里就不做一一介绍了。这里讲如何使用BWA MEM将质控合格的数据比对到参考基因组上。         BWA是一款基于BWT的快速比对工具，其由三个算法组成。这三个算法分别是：BWA backtrack, BWA SW and BWA MEM。其中，BWA MEM是最新的，其更快更准确，更适合用于人重数据分析。对于上述三种算法，首先需要使用索
生信：2：sam格式文件解读
热门推荐
genome_denovo的博客
12-04 5万+
第二章：生物信息分析第一节：解读sam格式文件1，SAM文件格式介绍SAM（The Sequence Alignment / Map format）格式，即序列比对文件的格式，详细介绍文档：http://samtools.github.io/hts-specs/SAMv1.pdfSAM文件由两部分组成，头部区和主体区，都以tab分列。 头部区：以’@’开始，体现了比对的一些总体信息。比如比对的SA
sam文件转换成bam文件
06-11
将SAM文件转换成BAM文件需要使用samtools工具，具体步骤如下： 1. 安装samtools工具，可以通过下面的命令进行安装： ``` sudo apt-get install samtools ``` 2. 使用samtools将SAM文件转换成BAM文件，命令如下： ``` samtools view -bS input.sam > output.bam ``` 其中，`input.sam`是输入的SAM文件名，`output.bam`是输出的BAM文件名。 3. 如果SAM文件中包含头部信息，还需要使用samtools将头部信息提取出来，命令如下： ``` samtools view -H input.sam > header.sam ``` 其中，`header.sam`是保存头部信息的SAM文件名。 4. 将头部信息和BAM文件合并成完整的BAM文件，命令如下： ``` samtools reheader header.sam input.bam > output.bam ``` 其中，`input.bam`是输入的BAM文件名，`output.bam`是输出的完整的BAM文件名。 这样，就可以将SAM文件转换成BAM文件了。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

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

余额充值