宏基因组测序流程(不完全版)

article/2025/11/7 1:01:54

文章目录

        • 所做工作
        • 收获
        • 宏基因组分析流程
          • Step1.去除宿主污染
          • Step2.去除接头序列
          • Step3.对序列进行进一步质控
          • Step4.对read进行进一步拼接(contig)
          • Step5.对contig进行orf预测
          • Step6.查找orf区翻译出来的氨基酸序列对应的蛋白质家族

所做工作

Step1.去除宿主污染

Step2.去除接头序列

Step3.对序列进⾏进⼀步质控

Step4.对read进⾏进⼀步拼接(contig)

Step5.对contig进⾏orf预测

Step6.查找orf区翻译出来的氨基酸序列对应的蛋⽩质家族

收获

  • 基因组学基本知识
  • 使用超算系统和在Linux下操作
  • 下载和使用软件
  • 理解各类常见数据格式的具体含义
  • 学会使用在线网站

宏基因组分析流程

概述:宏基因组也称微生物环境基因组,定义为生境中全部微小生物遗传物质的总和。目前主要指环境样品中的细菌和真菌的基因组总和。宏基因组学就是一种以环境样品中的微生群物体基因组为研究对象, 以功能基因筛选和/或测序分析为研究手段, 以微生物多样性、 种群结构、 进化关系、 功能活性、 相互协作关系及与环境之间的关系为研究目的的新的微生物研究方法。

其具体流程大致如下图所示:

在这里插入图片描述

Step1.去除宿主污染

BWA是一款将DNA序列mapping到参考基因组上的软件,例如比对到人类基因组。其由三个算法组成:BWA-backtrack,BWA-SW和BWA-MEM。

使用BWA整个比对过程主要分为两步,第一步建立参考基因组也就是宿主基因组的索引,第二步使用BWA MEM进行比对。

  • BWA建索引

    建索引有三种算法,通过参数-a is 、-a div和-a bwtsw进行选择。其中-a bwtsw对于短的参考序列是不工作的,必须要大于等于10Mb;-a is(效果和-a div是一样的)是默认参数,这个参数不适用于大的参考序列,必须要小于等于2G。

    • 格式:

      bwa index [-p prefix] [-a algoType] <in.db.fasta>
      
    • 实例:

      # 第一步,下载hg19基因组作为参考基因组
      wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
      tar zvfx chromFa.tar.gz
      cat *.fa > hg19.fa
      # 第二步,建立索引
      ../tools/bwa-0.7.17/bwa index -a bwtsw -p /mnt/disk1/yang/yanglab/wangjn/data/index/hg19 hg19.fa
      

在这里插入图片描述

  • 使用BWA MEM进行比对

    • 格式:

      bwa mem [-aCHMpP] [-t nThreads] [-k minSeedLen] [-w bandWidth] [-d zDropoff] [-r seedSplitRatio] [-c maxOcc] [-A matchScore] [-B mmPenalty] [-O gapOpenPen] [-E gapExtPen] [-L clipPen] [-U unpairPen] [-R RGline] [-v verboseLevel] db.prefix reads.fq [mates.fq]
      
      • 命令:
      ../../tools/bwa-0.7.17/bwa mem -t 16 -M  /mnt/disk1/yang/yanglab/wangjn/data/index/hg19 liweiling2016-10-25-shu-qian-2.read1_Clean.fastq.gz liweiling2016-10-25-shu-qian-2.read2_Clean.fastq.gz 1>liweiling2016-10-25-shu-qian-2.sam
      

在这里插入图片描述

  • 查看sam文件中和人类基因组比中的序列。
J00103:120:HFFJLBBXX:8:1101:6340:1613	163	chr11	55386856	0	34S19M97S	=	55386856	19	ATAAATCCGTTGCATGTACATTTTTGCTGTTATGGGTGTTCGTTTTTTTGTCCGCTGACTTGTTTGCCGGTTCATTCCGCACCTTGGGTATCAAGGAGGGATTGGGCAGTCGGCAAGTATTCCAGATAAACAAGGATTCAGCCGGATTTG	AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJFJJJJJJJJJJJJJJ<JJJJJFFJFJJJJJJJJJJFJJJJJJJJFF7FJJJJJJJJJJJJJJFFJJJ	NM:i:0	MD:Z:19	MC:Z:28S19M103S	AS:i:19	XS:i:0

比如说对于上面的这条结果,就是名称为J00103:120:HFFJLBBXX:8:1101:6340:1613的read与11号染色体的55386856位置比中。

  • 使用samtools去除宿主序列

    • 将sam文件转换为bam文件

      • 参数说明

        -b:以bam格式输出

        -u:以未压缩的bam格式输出,一般与linux命令配合使用

        -h:在sam输出中包含header信息

        -H:只输出header信息

        -f [INT]:只输出在比对flag中包含该整数的序列信息

        -F [INT]:跳过比对flag中含有该整数的序列

        -o [file]:标准输出结果文件

        -S : 默认下输入是 BAM 文件,若是输入是 SAM 文件,则最好加该参数,否则有时候会报错

      • 命令

        samtools view -bS liweiling2016-10-25-shu-qian-2.sam > liweiling2016-10-25-shu-qian-2.bam
        
    • 筛选出没有比中的序列

      • bam的flag图示,上面的标签都是2的n次方,这样的数有一个特点,就是随机挑选其中的几个,它们的和是唯一的。比如要筛选未比中的序列,由于4和8都代表未比中,因此筛选目标为4+8=12。

        在这里插入图片描述

        1 : 代表这个序列采用的是PE双端测序

        2: 代表这个序列和参考序列完全匹配,没有插入缺失

        4: 代表这个序列没有mapping到参考序列上

        8: 代表这个序列的另一端序列没有比对到参考序列上,比如这条序列是R1,它对应的R2端序列没有比对到参考序列上

        16:代表这个序列比对到参考序列的负链上

        32 :代表这个序列对应的另一端序列比对到参考序列的负链上

        64 : 代表这个序列是R1端序列, read1;

        128 : 代表这个序列是R2端序列,read2;

        256: 代表这个序列不是主要的比对,一条序列可能比对到参考序列的多个位置,只有一个是首要的比对位置,其他都是次要的

        512: 代表这个序列在QC时失败了,被过滤不掉了(# 这个标签不常用)

        1024: 代表这个序列是PCR重复序列(#这个标签不常用)

        2048: 代表这个序列是补充的比对(#这个标签具体什么意思,没搞清楚,但是不常用)

        • 命令
        samtools view -b -f 12 -F 256 liweiling2016-10-25-shu-qian-2.bam > liweiling2016-10-25-shu-qian-2-Unmapped.bam
        
        • 含义

        筛选出两端都没有比对上的序列,遇到不是主要比对的序列直接跳过。

    • 对bam进行排序

      • 说明

        samtools sort -n -m <in.bam> <out.prefix>

        -m 参数默认下是 500,000,000 即500M(不支持K,M,G等缩写)。处理大数据时,如果内存够用,则设置大点的值,以节约时间。-n 设定排序方式按short reads的ID排序。默认下是按序列在fasta文件中的顺序(即header)和序列从左往右的位点排序

      • 使用

        samtools sort -n liweiling2016-10-25-shu-qian-2-Unmapped.bam -o liweiling2016-10-25-shu-qian-2-Unmapped_sorted.bam
        
    • 使用bedtools将bam转换成fastq

    bedtools bamtofastq -i liweiling2016-10-25-shu-qian-2-Unmapped_sorted.bam -fq liweiling2016-10-25-shu-qian-2_removed_r1.fastq -fq2 liweiling2016-10-25-shu-qian-2_removed_r2.fastq
    

在这里插入图片描述

Step2.去除接头序列

可以采用SeqPrep程序合并重叠成单个较长读数的成对末端Illumina读数。也可以仅对adapter进行修剪,而不进行任何配对的末端重叠。

  • 必要的输入参数:
    • -f first read input fastqfilename #正链输入文件
    • -r second read input fastq filename #反链输入文件
    • -1 first read output fastq filename #正链输出文件
    • -2 second read output fastq filename #反链输出文件
  • 可选的输入参数:
    • -3 first read discarded fastq filename #正链丢弃的文件
    • -4 second read discarded fastq filename #反链丢弃的文件
    • -6#输入序列的格式,默认为phred+33,设置后为phred+64
    • -q Quality score cutoff for mismatches to be counted in overlap;default = 13 #质量分数
    • -L Minimum length of a trimmed or merged read to print it; default= 30 #最小长度
    • -O,-M,-N,-Q…
  • 关于adapter的控制参数:
    • -A #正链序列的adapter,默认值为AGATCGGAAGAGCGGTTCAG
    • -B #反链序列的adapter,默认值为AGATCGGAAGAGCGTCGTGT
  • 关于输出结果的控制参数:
    • -y #最大化质量分数
    • -s #指定将merge后的结果合并到哪个文件
    • -E,-x, -o, -m, -n指的是对输出结果进行进一步的细化操纵
../../tools/SeqPrep/SeqPrep -f liweiling2016-10-25-shu-qian-2_removed_r1.fastq -r liweiling2016-10-25-shu-qian-2_removed_r2.fastq -1 liweiling2016-10-25-shu-qian-2_removed_noadapter_r1.fq.gz -2 liweiling2016-10-25-shu-qian-2_removed_noadapter_r2.fq.gz -s liweiling2016-10-25-shu-qian-2_removed_noadapter.fq.gz
# 解压文件
gunzip liweiling2016-10-25-shu-qian-2_removed_noadapter.fq.gz
mv liweiling2016-10-25-shu-qian-2_removed_noadapter.fq liweiling2016-10-25-shu-qian-2_removed_noadapter.fastq

在这里插入图片描述

Step3.对序列进行进一步质控

FASTX-Toolkit是用于短读FASTA / FASTQ文件预处理的命令行工具的集合。FASTQ Information可用来执行图表质量统计和核苷酸分布,FASTQ/A Collapser可用来将FASTQ / A文件中的相同序列折叠成单个序列(同时保持读取计数)。FASTQ/A Clipper用来删除测序适配器/连接器。FASTQ Quality Filter用来根据质量过滤序列。要注意的一点是:不支持压缩格式的输入文件,不允许序列中存在N碱基,这样的序列会自动去除

  • fastx_clipper去除adapter
  • fastx_collapser比对基因组,去除重复序列
  • fastx_quality_stats对质量值进行统计
  • fastq_quality_boxplot_graph.sh绘制碱基质量分布盒式图
  • fastx_nucleotide_distribution_graph.sh绘制碱基分布图
fastx_collapser -v -i liweiling2016-10-25-shu-qian-2_removed_noadapter.fastq -o liweiling2016-10-25-shu-qian-2_removed_noadapter_collapsed.fastq && fastx_quality_stats -i liweiling2016-10-25-shu-qian-2_removed_noadapter_collapsed.fastq -o liweiling2016-10-25-shu-qian-2_removed_noadapter_collapsed.txt && fastq_quality_boxplot_graph.sh -i liweiling2016-10-25-shu-qian-2_removed_noadapter_collapsed.txt -o lliweiling2016-10-25-shu-qian-2_removed_noadapter_collapsed_quality.png -t "liweiling2016-10-25-shu-qian-2_removed_noadapter_collapsed" && fastx_nucleotide_distribution_graph.sh -i liweiling2016-10-25-shu-qian-2_removed_noadapter_collapsed.txt -o liweiling2016-10-25-shu-qian-2_removed_noadapter_collapsed_nuc.png -t "liweiling2016-10-25-shu-qian-2_removed_noadapter_collapsed"

查看图片

eog lliweiling2016-10-25-shu-qian-2_removed_noadapter_collapsed_quality.png
eog liweiling2016-10-25-shu-qian-2_removed_noadapter_collapsed_nuc.png

在这里插入图片描述
在这里插入图片描述

Step4.对read进行进一步拼接(contig)

细菌基因组组装工具之SPAdes SPAdes是由俄罗斯科学院圣彼得堡理工大学计算生物学实验室开发的基因组拼接工具。主要用于组装细菌基因组、二代(宏)基因组、(宏)转录组测序的拼接,也可用于一、二、三代测序的混合组装。是目前评价最好的拼接工具之一。

参数说明

  • Basic options:

    • -o 输出文件夹

    • –sc this flag is required for MDA (single-cell) data

    • –meta 这个是用于组装宏基因组**

    • –rna this flag is required for RNA-Seq data

    • –plasmid runs plasmidSPAdes pipeline for plasmid detection

    • –iontorrent this flag is required for IonTorrent data

    • –test runs SPAdes on toy dataset

    • -h/–help prints this usage message

    • -v/–version prints version

  • Input data:

    • –12 正向和反向合在一起的paired-end reads

    • -1 paired-end reads的正向序列

    • -2 paired-end reads的反向序列

    • -s unpaired reads的序列

  • Pipeline options:

    • –only-error-correction runs only read error correction (without assembling)

    • –only-assembler runs only assembling (without read error correction)

    • –careful tries to reduce number of mismatches and short indels

    • –continue continue run from the last available check-point

    • –restart-from restart run with updated options and from the specified check-point (‘ec’, ‘as’, ‘k’, ‘mc’)

    • –disable-gzip-output forces error correction not to compress the corrected reads

    • –disable-rr disables repeat resolution stage of assembling

  • Advanced options:

    • –dataset file with dataset description in YAML format
    • -t/–threads 线程数,默认是16个
    • -m/–memory RAM limit for SPAdes in Gb (terminates if exceeded) [default: 250]
    • –tmp-dir directory for temporary files [default: /tmp]
    • -k k-mer的设置数,如果是多个,中间要以逗号隔开,必须为奇数且不能大于128,即最大设置数 为127 [default: ‘auto’]
    • –cov-cutoff coverage cutoff value (a positive float number, or ‘auto’, or ‘off’) [default: ‘off’]

命令

/mnt/disk1/yang/yanglab/wangjn/tools/SPAdes-3.13.0-Linux/bin/spades.py -1 liweiling2016-10-25-shu-qian-2_removed_noadapter_r1.fq.gz -2 liweiling2016-10-25-shu-qian-2_removed_noadapter_r2.fq.gz -t 20 -o liweiling2016-10-25-shu-qian-2_removed_noadapter

在这里插入图片描述

Step5.对contig进行orf预测

下载了ORFfinder之后发现服务器下不能直接运行可执行文件,所以以线上预测代替。

网址:https://www.ncbi.nlm.nih.gov/orffinder/

  • 以Kmer值为21时计算得来的contig结果进行预测为例

原始序列

>NODE_6_length_670_cov_2.385208
TCCCACAAGATTGATATAATAATTTTCCAGACTTTCATCACATTCCTGCATAGCTTTAAT
CTTGCAATTTTTCTTGTCAAGAGTCAGCGCCAAATCTGTTACAGAAATACTGCTGAAAAT
ATCAGCCTCTGTTTGAGATAAAACAGTATACTCTACTCCCATTTCATCTAATACAAAGCT
CAAAGTTTTCGTGTCAGACACTTCTATGTGGATACATTTACGGCAGGCTTTTTCCAATTC
AACTGCACTGATTTCCTTTACAATACATCCTTTATCAATAAATCCAAAATGTGTTGCAAT
TTTTGACAACTCATCTAAAATGTGACTGGAAATCAAAACAGTAATCCCCTGCTCCCGATT
AAGCTTTAAAATCAATTCTCGAATTTCAATCATTCCTTGAGGGTCTAATCCATTTACCGG
TTCGTCTAAAATCAAAAAGTCCGGATTTCCTGCAAGTGCAACGGCAATTCCCAGTCTCTG
CCGCATACCCAGAGAAAAATTTTTTGCCTTTTTCTTACCGGTATTTTCCAGTCCCACAAG
CTTCAAAATATCAGTAATTCCTTCAAAAGACGGCATTCCTAGATACCGATATTGTGCTTT
CATATTATCTTCTGCCGTCATATCCAGATAAATAGCCGGTGTTTCCACCACAGCCCCTAT
TCTTTTCCTT

计算结果

在这里插入图片描述

Step6.查找orf区翻译出来的氨基酸序列对应的蛋白质家族

网址:http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi

上文计算出来的ORF2对应翻译的氨基酸序列

>lcl|ORF2
MTAEDNMKAQYRYLGMPSFEGITDILKLVGLENTGKKKAKNFSLGMRQRL
GIAVALAGNPDFLILDEPVNGLDPQGMIEIRELILKLNREQGITVLISSH
ILDELSKIATHFGFIDKGCIVKEISAVELEKACRKCIHIEVSDTKTLSFV
LDEMGVEYTVLSQTEADIFSSISVTDLALTLDKKNCKIKAMQECDESLEN
YYINLVG

查找结果

在这里插入图片描述

可以看到该段氨基酸序列由较大概率隶属于ABC_ATPase家族


http://chatgpt.dhexx.cn/article/pw1F4LVq.shtml

相关文章

一文读懂宏基因组分析套路

很多亲人感觉宏基因组的分析结果内容种类太多&#xff0c;根本学不过来。其实本质上并不复杂&#xff0c;只分为两类&#xff1a;物种组成和功能组成两大类&#xff0c;这是核心结果&#xff1b;再加上开头系统描述和结尾的讨论比较。通常会出现固定套路的4部分结构。 今天就从…

宏基因组测序

微基生物提供宏基因组测序分析服务。 宏基因组是指特定环境中全部微生物遗传物质的总和。宏基因组测序以特定环境中的整个微生物群落作为研究的对象&#xff0c;不需对微生物进行分离培养&#xff0c;而是提取环境微生物总DNA进行研究。其摆脱了传统研究中微生物分离培养的技术…

微生物菌群宏基因组研究技术分享

近年的研究热点集中于环境和生物体相互作用的微生物群体&#xff0c;而大量复杂的微生物群体存在培养困难&#xff0c;构成复杂&#xff08;包括细菌、古菌、真菌、原生生物、病毒甚至小型真核生物&#xff09;。因此如何用高通量精准的了解这些群体的构成&#xff0c;基因功能…

宏基因组分析-基于binning

一、介绍 宏基因组 ( Metagenome) 指特定环境下所有生物遗传物质的总和。它包含了可培养的和未可培养的微生物的基因。一般从环境样品中提取基因组DNA, 进行高通量测序&#xff0c;从而分析微生物多样性、种群结构、功能信息、与环境之间的关系等。 宏基因组的分析目前主要包…

宏基因组生信分析方法介绍

随着高通量技术的发展&#xff0c;宏基因组学&#xff08;metagenomics&#xff09;已经成为研究微生物群落物种及功能的前沿科学&#xff0c;在肠道微生物、环境微生物等研究领域具有广泛应用。宏基因组学通过对微生物群落全部DNA进行高通量测序&#xff0c;将测序序列与公共数…

宏基因组分析-基于组装

一、介绍 宏基因组 ( Metagenome) 指特定环境下所有生物遗传物质的总和。它包含了可培养的和未可培养的微生物的基因。一般从环境样品中提取基因组DNA, 进行高通量测序&#xff0c;从而分析微生物多样性、种群结构、功能信息、与环境之间的关系等。 宏基因组的分析目前主要包…

逻辑回归的常见问题

文章目录 逻辑回归概括逻辑回归的基本假设逻辑回归的损失函数交叉熵损失函数的原理交叉熵损失函数的直观理解logistic回归和线性回归的关系线性回归和逻辑回归的异同逻辑回归是线性模型吗类别的分界线是什么逻辑回归输出值是概率吗逻辑回归也可以处理多分类逻辑回归的求解方法梯…

SQL逻辑运算

SQL逻辑运算符 SQL 逻辑运算符逻辑运算符 ALL逻辑运算符 AND逻辑运算符 ANY逻辑运算符 BETWEEN逻辑运算符 EXISTS逻辑运算符 IN逻辑运算符 LIKE逻辑运算符 NOT逻辑运算符 OR逻辑运算符 IS NULL SQL 逻辑运算符 运算符描述ALL所有运算符用于比较的值到另一个值组中的所有值ANDA…

逻辑英语-写作

指日可待 In this way, a better tomorrow will not be a question of if, but when.1、shoulder 我们应当承担起保护环境的重任。 We must shoulder the liability of -----2. Word 我对你的感谢已经溢于言表 My thanks to you is beyond words.3. house 容纳 A wise man she…

逻辑英语公式R1+R2

一、主谓宾 1.1、主语 一般是名词&#xff0c;一般是动作的发出者&#xff0c;或者是被描述的对象 1.2、宾语 一般也是名词&#xff0c;一般是主语作用到的对象。 1.3、谓语 如果没有谓语那么这个世界就是静态&#xff0c;谓语就是来让这个主语和谓语产生关系。 常见的是…

串口USART和UART

串口通信&#xff1a; UART是通用串行数据总线&#xff0c;用于异步通信。该总线双向通信&#xff0c;可以实现全双工传输和接受。UART主要用于主机与辅助设备通信。 UART的功能计算器内部采用并行数据&#xff0c;不能直接把数据发到Modem&#xff0c;必须经过UART整理才能进…

【通信接口】UART、IIC、SPI

目录 一、预备知识 1、串行与并行 2、单工与双工 3、波特率 二、UART 三、IIC 四、SPI &#xff08;一对一、一对多&#xff09; 五、IIC、SPI异同点 参考文章&#xff1a;这些单片机接口&#xff0c;一定要熟悉&#xff1a;UART、I2C、SPI、TTL、RS232、RS422、RS485…

基于FPGA的UART接口设计

一、顶层设计思路&#xff1a; UART即通用异步收发传输接口&#xff08;Universal Asynchronous Receiver/Transmitter&#xff09;&#xff0c;简称串口&#xff0c;是一种常用的通信接口&#xff0c;其协议原理就不赘述了&#xff0c;不了解的可以自己查阅资料。&#xff08;…

单片机通信接口:UART、I2C、SPI、TTL、RS232、RS422、RS485、CAN、USB

参考资料&#xff1a; 这些单片机接口&#xff0c;一定要熟悉&#xff1a;UART、I2C、SPI、TTL、RS232、RS422、RS485、CAN、USB、SD卡 秒懂所有USB接口类型&#xff0c;USB接口大全 1. UART UART(通用异步收发器)指的是一种物理接口形式(硬件)。 UART是异步&#xff0c;全双…

通信接口:UART、I2C、SPI、TTL、RS232、RS422、RS485、CAN、USB

1. UART UART(通用异步收发器)指的是一种物理接口形式(硬件)。 UART是异步&#xff0c;全双工串口总线。它比同步串口复杂很多。有两根线&#xff0c;一根TXD用于发送&#xff0c;一根RXD用于接收。 UART的串行数据传输不需要使用时钟信号来同步传输&#xff0c;而是依赖于发送…

USB,串口(RS232、RS485),UART接口

USB转串口即实现计算机USB接口到通用串口之间的转换。为没有串口的计算机提供快速的通道&#xff0c;而且&#xff0c;使用USB转串口设备等于将传统的串口设备变成了即插即用的USB设备。作为应用最广泛的USB接口&#xff0c;每台电脑必不可少的通讯接口之一&#xff0c;它的最大…

Uart接口的详细解释

我面试的时候一般喜欢问应聘者一个问题&#xff1a;UART与RS232/RS485的区别与联系&#xff1f;很多人对于这个问题答得都不是很好。还有些人压根就没有想过这个问题&#xff0c;一直认为他们是同一个东西&#xff0c;就是咱们俗称的串口。 我刚入嵌入式的大门时&#xff0c;对…

UART接口说明

逼近年关事情多&#xff0c;少了更新。今天冒个泡。说下UART通信接口。 UART扫盲 前面做了SPI和I2C&#xff0c;前两者一个是摩托&#xff0c;一个是飞利浦背书&#xff0c;简单好理解。这个UART就相对复杂一点&#xff0c;全称universal Asynchronous Receiver/Transmitter …

UART接口详解

文章目录 简介硬件接线RS232RS485RS232和RS485比较 通信原理uart和usart的区别实例针对STM32的串口数据位特点&#xff0c;改成对应PC的串口数据校验当使用9600波特率的时候&#xff0c;通讯稳定&#xff0c;当使用115200波特率的时候&#xff0c;通讯变得不稳定。 简介 UART全…

UART接口介绍

0 Preface/Foreword UART是Universal Asynchronous Receiver and Transmitter简称&#xff0c;中文为 通用异步接收和发送器&#xff0c;是常用的串行通讯接口。 RS-232&#xff1a;RS-232标准接口&#xff08;aka. EIA RS-232&#xff09;是常用的串行通信接口标准之一&#…