SMURF-Science封面文章使用的16S新流程

最新推荐文章于 2023-07-14 14:06:45 发布
最新推荐文章于 2023-07-14 14:06:45 发布
阅读量662 收藏 3
点赞数 1
分类专栏： 生物信息 文章标签： SMURF 16S
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/zd200572/article/details/113039048
版权

在这里插入图片描述
肠道微生物是近两年的研究热点，但是去年登上Science封面的是一篇研究肿瘤中的微生物的文章，另人眼前一亮，有些肿瘤即使没有与外界环境相通，也是有微生物的存在的。外行看热闹，内行要看看他是具体怎么进行研究的。
在这里插入图片描述
首先是研究手段，并不是宏基因组，是16S，估计是由于肿瘤中的微生物含量过少，多数不能满足宏基因组的建库所需DNA的量。然后，作者是用了一种不同于常规16S的研究手段进行的，扩增并测序了5段V区（68%的长度），然后合并分析的，作者称之为SMURF的方法流程，认为这个方法是接近于三代16S全长的物种分辨率的，并给出了参考文献。可应用于任何一组扩增子，从而实现有效的下游分析。
在这里插入图片描述
鉴于三代测序的成本居高不下，这个方法还是有一定市场的，二代的白菜价格，获得三代的结果，何乐而不为呢？有这么好的学习资源，我们当然要学习一下嘛。建库实验方面并没有多大的问题，我们主要来看下数据分析的部分。
在这里插入图片描述
算法的运行方式有两种，matlab里面运行，类似R语言，或者依赖于MCR库，不需要安装matlab（类似于R语言的运行方式吧），我选择了后者，毕竟matlab收费的。
在这里插入图片描述
还有个超级可爱的图标，github上看不到，拖回国内才发现。

软件安装与配置

提前声明，这个脚本会报错，不能使用，如想使用可采用qiime2插件进行。以下内容可无视。


# 下载依赖的matlab MCR平台（作者使用matlab写的分析软件）虽然这是美帝的，但是学习先进技术嘛！
mkdir MCR && cd MCR
wget https://www.mathworks.com/supportfiles/downloads/R2014a/deployment_files/R2014a/installers/glnxa64/MCR_R2014a_glnxa64_installer.zip
unzip MCR_R2014a_glnxa64_installer.zip
# 安装下，必须有图形界面才能安装成功。
sudo ./install
# 设置环境变量，这里是临时的，所以退出终端后埼再添加一次
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/MATLAB/MATLAB_Compiler_Runtime/v83/runtime/glnxa64:/usr/local/MATLAB/MATLAB_Compiler_Runtime/v83/bin/glnxa64:/usr/local/MATLAB/MATLAB_Compiler_Runtime/v83/sys/os/glnxa64
export XAPPLRESDIR=/usr/local/MATLAB/MATLAB_Compiler_Runtime/v83/X11/app-defaults
# 报了个这个，因为找不到图形界面，算法不需要，应该不影响
#Exception in thread "main" java.lang.InternalError: Can't connect to X11 window server using ':0' as the value of the DISPLAY variable.
cd ..
# 下载脚本
# https://github.com/NoamShental/SMURF.git
# 我把它拉到了gitee上，克隆速度快上许多，特别是这个，因为挺大的
git clone https://github.com/NoamShental/SMURF.git
cd SMURF
# 数据库准备
cat ./Green_Genes_201305/unique_up_to_3_ambiguous_16S/GreenGenes_201305_unique_up_to_3_ambiguous_16S.fasta.gz*> ./Green_Genes_201305/unique_up_to_3_ambiguous_16S/Green_Genes_201305_unique_up_to_3_ambiguous_16S.fasta.gz
gunzip ./Green_Genes_201305/unique_up_to_3_ambiguous_16S/Green_Genes_201305_unique_up_to_3_ambiguous_16S.fasta.gz

运行与结果

运行好像一条命令就行了，前提是配置好引物等参数。
需要修改的参数：

% ********************** GENERAL PARAMETERS ********************
base_samples_dir = '/';
...
% ********************** SAMPLE PREP PARAMETERS ********************
% Set the 16S reference DB
uniS16_dir = './Green_Genes_201305/unique_up_to_3_ambiguous_16S';
db_filename = 'Green_Genes_201305_unique_up_to_3_ambiguous_16S';
# 其他参数，不确定是否需要
 vi Configs/db_params_script.m
 #把 ../ 替换为./或者在Standalone文件夹运行，不需要改
vi Configs/adhoc_db_params_script.m

运行啦

chmod +x ./StandaloneVersion/SMURF_lin
 time ./StandaloneVersion/SMURF_lin ./Configs/compiled_params_script.m

当然，示例文件肯定不会报错，很轻松出结果嘛。


```bash
time  ./StandaloneVersion/SMURF_lin ./Configs/compiled_params_script.m
Doing quality filters
Part 1/1 - Block 1/5
Part 1/1 - Block 2/5
Part 1/1 - Block 3/5
Part 1/1 - Block 4/5
Part 1/1 - Block 5/5
Number of reads: 472350
Percent of long enough reads: 0.94713
Percent of good reads: 0.91592
Counting fasta write: 1
Elapsed time is 9.831863 seconds.
Mapped to primers 82% of unique reads
Mapped to primers 97% of read counts
regions_files = 
6x1 struct array with fields:
    name
    date
    bytes
    isdir
    datenum
ans =
./Green_Genes_201305/unique_up_to_3_ambiguous_16S_amp6Regions_2mm_RL130/GreenGenes_201305_unique_up_to_3_ambiguous_16S_amp6Regions_2mm_RL130_region1.mat
Loading bacterial DB for region 1 out of 6 from original region 1
Loading bacterial DB for region 2 out of 6 from original region 2
Loading bacterial DB for region 3 out of 6 from original region 3
Loading bacterial DB for region 4 out of 6 from original region 4
Loading bacterial DB for region 5 out of 6 from original region 5
Loading bacterial DB for region 6 out of 6 from original region 6
Region 1 out of 6
Keep high freq: 28% of reads
Keep high freq: 91% of counts
Building matrix M
Building matrix Q
 --------------------------------------------
...
--------------------------------------------
Region 6 out of 6
Keep high freq: 2% of reads
Keep high freq: 89% of counts
Building matrix M
Building matrix Q
--------------------------------------------
Region 1 out of 6
Keeping reads matched to DB: 95% of reads
Keeping reads matched to DB: 98% of counts
--------------------------------------------
...
--------------------------------------------
Region 6 out of 6
Keeping reads matched to DB: 97% of reads
Keeping reads matched to DB: 100% of counts
--------------------------------------------
Filter out columns (bacteria)
Normalize frequency counts
Build matrix A_L2
Iter:4674. Error reduction of X (L1 norm): 9.7149e-07
Total iterations time: 60.4761
Error using main_multiple_regions (line 34)
Not enough input arguments.
Error in main_smurf (line 36)

MATLAB:minrhs

real	9m1.616s
user	11m32.462s
sys	0m38.527s

所以它是报错了，这个代码有问题，但是github上发现了另一个实现方式，qiime2 Sidle插件，也可以做到，切换工具，继续学习。

zd200572
关注
  • 1
    点赞
  • 踩
  • 3
    收藏
    觉得还不错? 一键收藏
  • 5
    评论
专栏目录
Surf算法特征点检测与匹配
牧野的博客
09-20 11万+
Speeded Up Robust Features（SURF，加速稳健特征），是一种稳健的局部特征点检测和描述算法。最初由Herbert Bay发表在2006年的欧洲计算机视觉国际会议（Europen Conference on Computer Vision，ECCV）上，并在2008年正式发表在Computer Vision and Image Understanding期刊上。 Surf
SMURF(5R)-Science封面文章使用的16S新流程(二)
zd200572的博客
01-30 690
前面介绍的SMURF流程的运行以失败告终了，不过这个是这篇文章的参考方法，至于这篇文章改进过的方法，还没有试过，这就试一下，顺便考虑是否能把6区的移植过来，搞个6R呢，可能，算法上有略微的区别，毕竟这篇Science研究的是肿瘤中的含量很少的微生物，用了严格的去污染策略，不管怎样，试试吧！ 1、环境准备 类似上次那个流程，更加简单了些，只需要安装解压下。 # 安装MCR,这次是新版本的9.7,重要的事说三遍，必须有图形界面gui，否则会安装失败 #必须有图形界面gui，必须有图形界面gui #　下载地址，速
【阅读论文】时间序列数据清洗:一项调查
zijin-jdd的博客
06-17 1799
误差在时间序列数据中普遍存在，在工业领域尤为普遍。错误的数据无法存储在数据库中，导致数据资产丢失。目前，为了处理这些包含错误的时间序列，除了保留原始错误数据、丢弃错误数据和手动检查错误数据外，我们还可以使用数据库中广泛使用的清洗算法对时间序列数据进行自动清洗。本调查提供了时间序列数据清洗技术的分类，并全面回顾了每种类型的最先进的方法。此外，我们还总结了来自研究和行业的数据清理工具、系统和评估标准。最后，我们强调了时间序列数据清洗的可能方向。关键词：数据清洗，数据质量，时间序列。
SMURF使用经验
最新发布
xiaozheng1213的博客
07-14 263
5R 16S比对
SMURF流程之q2-sidle（四）-- 序列重建
zd200572的博客
02-06 518
‎SMURF 算法的核心是基于基于 kmer 的短区域重建到全长框架中。有两个步骤，首先是ASV在单个区域基于kmer进行比对，然后完整的序列集组装成重建的计数表。 区域比对 第一步是每个区域把序列比对到数据库，使用 align-regional-kmers 命令，我们前面使用--kmer-db-fp选项设置了数据库，使用 --rep-seq-fp选项传递ASV序列，最后是区域定义，来自前面你给区域起的别名，要完全一致。比对是一个开心的可并行任务，我们可以通过多多线程提升性能(--p-n-workers参
Smurf攻击的原理
热门推荐
芝兰生于深谷，不以无人而不芳
05-07 2万+
原理：攻击者向网络广播地址发送ICMP包，并将回复地址设置成受害网络的广播地址，通过使用ICMP应答请求数据包来淹没受害主机的方式进行，最终导致该网络的所有主机都对次ICMP应答请求作出答复，导致网络阻塞。更加复杂的Smurf攻击攻击将源地址改为第三方受害者，最终导致第三方崩溃。 防范措施： 1 配置路由器禁止IP广播包进网 2 配置网络上所有计算机的操作系统，禁止对目标地址为广播地址
Fat Smurf-crx插件
04-06
"Fat Smurf-crx插件"是一个专门为网络浏览器设计的扩展程序，主要功能是将用户在浏览网页时遇到的所有图像替换为一个特定的卡通角色——胖爸爸蓝精灵。这个插件以英语为操作语言，旨在增加用户的浏览乐趣，特别是在...
Smurf SoundFont Editor-开源
05-09
**Smurf SoundFont Editor——一个开源的音频编辑工具** Smurf SoundFont Editor是一款专为音乐制作人和声音设计师设计的开源软件，它允许用户编辑和定制SoundFont文件。SoundFont是一种流行的音频格式，广泛用于...
matlab匹配滤波代码-5R:基于5R协议的SMURF测序结果的重建
05-24
使用我们的短多区域框架（SMURF）算法（Fuks等人2018），将来自这些区域的结果读取合并为一个一致的解决方案。 提供了两个选项：Matlab代码和Linux的独立可执行版本。 两种方法的输入都是相同的，因此下面将针对...
Smurf Combat:Smurf Combat BBS DOOR Game（大约 1993-1995）-开源
06-01
Smurf Combat BBS DOOR Game（大约 1993-1995） 版权所有 (C) 1993, 1994, 1995, 2019 Laurence Maar ... 您可以根据自由软件基金会发布的 GNU 通用公共许可证的条款重新分发和/或修改它； 许可证的第 2 版，或（由您...
SMURF和D.R.DOS
starshift的专栏
03-10 1788
一.SMURF攻击的概念及效果1.概念:Smurf攻击的命名，是从smurf这个攻击工具来的。Smurf是一种在局域网中的攻击手段，它的作用原理是基于广播地址与回应请求的，我们都知道，在局域网中，计算机被分配的ip地址的最后一位，是在1-254之间，那么最后那个255是留下来干什么的呢？其实带有这个255的地址，比如:192.168.0.255，就是这个子网的广播地址你可以这样试试看:运行->CM
smurf攻击
金溪的博客
10-23 8247
    ICMP ECHO请求包用来对网络进行诊断，当一台计算机接收到这样一个报文后，会向报文的源地址回 应一个ICMP ECHO REPLY。一般情况下，计算机是不检查该ECHO请求的源地址的，因此，如果一个恶意 的攻击者把ECHO的源地址设置为一个广播地址，这样计算机在恢复REPLY的时候，就会以广播地址为目的 地址，这样本地网络上所有的计算机都必须处理这些广播报文。如果攻击者发送的ECH...
smurf分布式攻击（DDOS）
一个帅气潇洒的豆子
01-23 4979
Smurf攻击是一种逼毒攻击，以最初发动这种攻击的程序“Smurf”来命名。这种攻击方法结合使用了IP欺骗和ICMP回复方法使大量网络传输充斥目标系统，引起目标系统拒绝为正常系统进行服务。(引自百度) The Smurf attack is a distributed denial-of-service attack in which large numbers of Internet Co
SMURF流程之q2-sidle（－）及Qiime2安装加速
zd200572的博客
01-24 510
前面说到Science封面文章用的16S数据分析流程有qiime2的插件版本，可以解决基于matlab MCR standalone版本的报错，于是实践一下！https://github.com/jwdebelius/q2-sidle conda的安装就不表了，教程挺多的。 环境准备 安装qiime2-2020.11 作者说只测试了兼容这个版本，于是就装这个啦！ # 激活环境 source ~/data_home/Miniconda3/bin/activate # 下载配置文件 wget https://d
SMURF流程之q2-sidle（三）--reads准备
zd200572的博客
02-06 691
完成了前面的数据库准备，下面就是reads的准备，基本过程就是把reads拆成对应不同引物的几个部分，后面再重建合并在一起啦。首先声明，这个方法还在开发和完善之路，最近一次更新在这个月，可能结果会有变动，应该说还处于beta版本中，不建议在生产环境中使用。 这里就有几种情况啦，一种是已经每个样本每个V区拆好的数据，另一种是每个样本几个V区混在一起的数据，或者完全没拆的数据。这里根据SMURF的示例，按第二种情况进行，应该是最常见的情况。下面是具体步骤： Reads准备 尽管SMURF依赖于质控过滤，还是推荐
pytorch实现LSTM训练模型，使用NSS-KDD数据集
04-04
首先，我们需要加载NSS-KDD数据集。可以从以下链接下载数据集：http://www.unb.ca/cic/datasets/nsl.html 在代码中，我们使用pandas库来读取csv文件，并将标签转换为数字。我们还需要使用sklearn库中的train_test_split函数来将数据集分为训练集和测试集。 ```python import pandas as pd from sklearn.model_selection import train_test_split # 加载数据集 df = pd.read_csv('KDDTrain+.csv') # 将标签转换为数字 df['label'] = df['label'].map({'normal': 0, 'neptune': 1, 'warezclient': 2, 'ipsweep': 3, 'portsweep': 4, 'teardrop': 5, 'nmap': 6, 'satan': 7, 'smurf': 8, 'pod': 9, 'back': 10, 'guess_passwd': 11, 'ftp_write': 12, 'multihop': 13, 'rootkit': 14, 'buffer_overflow': 15, 'imap': 16, 'warezmaster': 17, 'phf': 18, 'land': 19, 'loadmodule': 20, 'spy': 21, 'perl': 22, 'saint': 23, 'mscan': 24, 'apache2': 25, 'snmpgetattack': 26, 'processtable': 27, 'httptunnel': 28, 'ps': 29, 'snmpguess': 30, 'mailbomb': 31, 'named': 32, 'sendmail': 33, 'xterm': 34, 'worm': 35, 'xlock': 36, 'xsnoop': 37, 'sqlattack': 38}) # 分割数据集为训练集和测试集 X_train, X_test, y_train, y_test = train_test_split(df.drop('label', axis=1), df['label'], test_size=0.2) ``` 接下来，我们需要对数据进行预处理。我们使用sklearn的StandardScaler函数对数据进行标准化。此外，我们还需要将数据转换为PyTorch张量。 ```python import torch from sklearn.preprocessing import StandardScaler # 对数据进行标准化 scaler = StandardScaler() X_train = scaler.fit_transform(X_train) X_test = scaler.transform(X_test) # 将数据转换为PyTorch张量 X_train = torch.tensor(X_train).float() X_test = torch.tensor(X_test).float() y_train = torch.tensor(y_train.values) y_test = torch.tensor(y_test.values) ``` 现在，我们可以构建LSTM模型。在这个例子中，我们使用两层LSTM和一个全连接层。我们还需要定义一个损失函数和一个优化器。 ```python import torch.nn as nn import torch.optim as optim # 定义LSTM模型 class LSTM(nn.Module): def __init__(self, input_size, hidden_size, num_layers, num_classes): super(LSTM, self).__init__() self.hidden_size = hidden_size self.num_layers = num_layers self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True) self.fc = nn.Linear(hidden_size, num_classes) def forward(self, x): h0 = torch.zeros(self.num_layers, x.size(0), self.hidden_size).to(device) c0 = torch.zeros(self.num_layers, x.size(0), self.hidden_size).to(device) out, _ = self.lstm(x, (h0, c0)) out = self.fc(out[:, -1, :]) return out # 定义模型参数 input_size = X_train.shape[1] hidden_size = 128 num_layers = 2 num_classes = 39 learning_rate = 0.001 # 定义模型、损失函数和优化器 model = LSTM(input_size, hidden_size, num_layers, num_classes).to(device) criterion = nn.CrossEntropyLoss() optimizer = optim.Adam(model.parameters(), lr=learning_rate) ``` 现在，我们可以开始训练模型。我们需要将数据分批次传递给模型，并在每个批次之后更新模型的权重。 ```python # 训练模型 num_epochs = 10 batch_size = 64 for epoch in range(num_epochs): for i in range(0, len(X_train), batch_size): inputs = X_train[i:i+batch_size].to(device) targets = y_train[i:i+batch_size].to(device) # 前向传播 outputs = model(inputs) # 计算损失 loss = criterion(outputs, targets) # 反向传播和优化 optimizer.zero_grad() loss.backward() optimizer.step() # 每个epoch打印损失 print(f'Epoch {epoch+1}/{num_epochs}, Loss: {loss.item():.4f}') ``` 最后，我们可以使用测试集评估模型的性能。 ```python # 评估模型 with torch.no_grad(): correct = 0 total = 0 for i in range(0, len(X_test), batch_size): inputs = X_test[i:i+batch_size].to(device) targets = y_test[i:i+batch_size].to(device) # 前向传播 outputs = model(inputs) # 预测类别 _, predicted = torch.max(outputs.data, 1) total += targets.size(0) correct += (predicted == targets).sum().item() print(f'Test Accuracy: {100 * correct / total:.2f}%') ``` 完整代码如下：
评论 5
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值