# 1.Mapping * 二代测序数据种类繁多，分析方法也各有差异。在待分析的物种已有参考基因组的情况下，传统的高通量数据分析流程中通常都会把测序数据mapping回参考基因组(对RNA-seq在有一些情况下有的流程也会考虑向转录组mapping)，再用mapping的结果（通常是一个bam格式的文件）进行后续的分析。 ## 0)准备运行环境 进入我们之前使用的docker容器: ```bash docker exec -it bioinfo_tsinghua /bin/bash # 进入本节课的工作目录 cd /home/test/mapping ``` ## 1) 流程 ![](/files/-LPVvEulZ6EGwv9YRMFw) ## 2) 文件格式 ### 2a) fastq文件 * fastq文件是存储二代测序测出的reads序列的最常见的一种文件格式。 * fastq文件除了read的序列信息之外，还记录了每一个碱基的测序质量信息 * fastq文件中，每四行对应一个read，下面提供了一个简单的例子。 * 第1行以"@"开头，"@"后面记录了read id，read id后面还可以空一格，再加上一些相关的描述信息 * 第2行记录着read序列 * 第3行以"+"开头，通常包含一些相关的描述信息(如果为空，也需要一个"+"作为占位符) * 第4行为ASCII码表示的每个碱基的测序质量，长度和第2行相等 ``` @EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC + ;;3;;;;;;;;;;;;7;;;;;;;88 ``` * fastq文件可以很容易的转换成fasta文件，但这一步会丢失测序质量的信息。大部分mapping的软件也支持fasta文件格式作为输入 * 对于双端测序，两个end的reads通常放在两个fastq文件中，在这两个fastq文件中，对应同一个read pair的两个read次序是一致的(例如对于第一个fastq文件中的第43个record对应的read，与它paired的read位于第二个fastq文件的第43个record)。 ### 2b) sam文件 * sam/bam 文件是存储二代测序数据mapping结果的最常见的文件格式, mapping的过程可以简单的理解为一个从fastq文件产生bam文件的过程 * sam文件为纯文本文件，bam文件为压缩后的二进制文件。在实际工作中bam文件更常用(因为比较节约存储) * sam/bam文件一般由两部分组成: 第一部分称为header，它包含了reference序列的名称及长度，mapping结果的排序方式，mapping使用的软件等信息；第二部分为sam/bam文件的主体，包含了reads mapping的信息 * 在不考虑multiple alignment(一个read被比对到多个位置)和chimeric alignment(一个read的不同区域被比对到了不同的染色体等)的情况下，一条read mapping的结果通常对应着sam/bam文件主体部分的一行 * sam/bam文件可以包括unmapped reads，也可以不包括unmapped reads * pair end reads mapping的结果通常存储在一个bam文件中 * sam/bam 文件的具体定义相对比较复杂，对细节感兴趣的同学可以参考 * samtools为我们操纵sam/bam文件提供了很多非常方便的功能，我们将在samtools/bedtools一节进行具体介绍 ## 3) unspliced reads mapping * 在向基因组mapping各个物种的DNA-seq(以及原核生物的RNA-seq)数据，或向转录组mapping各个物种的RNA-seq数据时，一般不需要考虑RNA splicing的问题。 * 对于这样的任务，[bowtie](http://bowtie-bio.sourceforge.net/index.shtml),[bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml)和[bwa](http://bio-bwa.sourceforge.net/)是比较常用的工具。 * 我们这里主要介绍[bowtie](http://bowtie-bio.sourceforge.net/index.shtml)的基本使用。 ```bash # 将fasta文件中THA1.fa的reads mapping到酵母基因组上 bowtie -v 2 -m 10 --best --strata BowtieIndex/YeastGenome -f THA1.fa -S THA1.sam # -v 2: 最多容许2个mismatch # -m 10: 只输出可以map到不超过10个位置的reads mapping的结果 # --best --strata: 只汇报最好的一个hit,两个参数需要同时指定 # BowtieIndex/YeastGenome: 酵母的bowtie index,可以从https://bowtie-bio.sourceforge.net/manual.shtml下载，也可以用bowtie-build从基因组文件自己建立 # -f THA1.fa: 输入为fasta文件，路径为THA1.fa # -S THA1.sam: 输出文件名为THA1.sam，格式为sam文件 # 将fastq文件e_coli_1000_1.fq中的reads mapping到大肠杆菌基因组上 bowtie -v 1 -m 10 --best --strata bowtie-src/indexes/e_coli -q e_coli_1000_1.fq -S e_coli_1000_1.sam # -v 1: 最多容许1个mismatch # -m 10: 只输出可以map到不超过10个位置的reads mapping的结果 # --best --strata: 同上 # bowtie-src/indexes/e_coli 大肠杆菌的bowtie index,可以从https://bowtie-bio.sourceforge.net/manual.shtml下载，也可以用bowtie-build从基因组文件自己建立 # -q e_coli_1000_1.fq: 输入为fastq文件 # -S e_coli_1000_1.sam: 同上 ``` {% hint style="info" %} * bowtie既可以用来mapping单端测序数据，也可以用来mapping双端测序数据。如果希望了解mapping双端测序的命令，请参考bowtie的文档: * bowtie在mapping中是不考虑insertion和deletion的，而bowtie2和bwa会考虑这些情况。有兴趣了解的同学请参考它们的文档，了解它们的使用方法。 {% endhint %} ## 4) Genome Browser see [1.1-genome-browser](/teaching/part-iii.-ngs-data-analyses/1.mapping/1.1-genome-browser.md) ## 5) 延伸阅读 ### 5a) 基因组学的常用文件格式 详见 。 ### 5b) spliced mapping * 对于真核生物的RNA-seq数据的处理，我们通常需要考虑RNA splicing的问题。 * 有很多软件是专门针对向真核基因组mapping RNA-seq数据设计的，[hisat2](http://daehwankimlab.github.io/hisat2/)和[STAR](https://github.com/alexdobin/STAR)目前来说是两个非常常用的软件。[tophat](https://ccb.jhu.edu/software/tophat/index.shtml)曾经也是一个流行的工具，但现在用的人已经不多了。 * mapping RNA-seq数据的工具通常还需要提供一个gtf注释文件以获得splicing位点的信息，再向index mapping测序数据，产生bam文件。 * 这里我们也提供了利用[STAR](https://github.com/alexdobin/STAR)软件向拟南芥基因组mapping RNA-seq reads的例子。 * 在docker镜像中下载STAR的release并解压(在宿主机下载后，通过共享目录在docker中使用当然也是可以的): ```bash wget https://github.com/alexdobin/STAR/archive/refs/tags/2.7.10a.zip unzip 2.7.10a.zip ``` * 下载所需的数据[STAR.mapping.tar.gz](https://cloud.tsinghua.edu.cn/d/429647f35add41338d1a/)放入容器内的当前工作目录并解压: ```bash tar xvzf STAR.mapping.tar.gz ``` {% hint style="info" %} * 为了降低计算开销，供大家学习之用，我们这里使用的数据是从一个拟南芥RNA-seq数据mapping到叶绿体的Reads中sampling出来的，参考基因组和基因组注释也只包含了拟南芥叶绿体的序列。 * 这里提供的是双端测序的数据，单端测序数据的mapping非常类似，请自行参看STAR的文档。大家可以用这里双端测序的一个end(如`ath_1.fastq`)当作单端测序数据进行尝试。 {% endhint %} * 建立STAR index: ```bash mkdir tair10.Pt.STARindex STAR-2.7.10a/bin/Linux_x86_64_static/STAR --runMode genomeGenerate --genomeFastaFiles data/tair10.Pt.fa --sjdbGTFfile data/tair10.Pt.gtf --genomeDir tair10.Pt.STARindex --genomeSAindexNbases 7 # genomeSAindexNbases 是STAR suffix array pre-indexing 的字符串长度, 默认为14 # 对于比较小的基因组(例如我们这里只用到了叶绿体序列，可以理解成一个很小的基因组)一般会设成min(14, log2(GenomeLength)/2 - 1) ``` * mapping ```bash mkdir output STAR-2.7.10a/bin/Linux_x86_64_static/STAR --runMode alignReads --genomeDir tair10.Pt.STARindex --readFilesIn data/ath_1.fastq data/ath_2.fastq --outFileNamePrefix output/ath.aligned ``` * 在输出目录中，`output/ath.alignedAligned.out.sam`是mapping产生的sam文件，`output/ath.alignedSJ.out.tab`是splice junction的信息，其余两个是STAR输出的日志文件，包括mapping结果的统计信息等。 {% hint style="info" %} * 我们这里只介绍了如何用STAR的默认参数进行reads mapping，没有涉及如何结合实际需求调整各种参数。STAR是一个非常灵活的工具，它提供了大量的参数和多种功能，有兴趣的同学请参考文档 * 我们这里只介绍了bowtie和STAR，也鼓励大家用我们这里提供的数据，用不同的aligner(bowtie2,bwa,hisat2等)进行各种尝试 {% endhint %} ## 6) Homework * （1）请阐述bowtie中利用了 BWT 的什么性质提高了运算速度？并通过哪些策略优化了对内存的需求？ * （2）用bowtie将 `THA2.fa` mapping 到 `BowtieIndex/YeastGenome` 上，得到 `THA2.sam`，统计mapping到不同染色体上的reads数量(即统计每条染色体都map上了多少条reads)。 * （3）查阅资料，回答以下问题: * （3.1）什么是sam/bam文件中的"CIGAR string"? 它包含了什么信息? * （3.2）"soft clip"的含义是什么，在CIGAR string中如何表示？ * （3.3）什么是reads的mapping quality? 它反映了什么样的信息? * （3.4）仅根据sam/bam文件的信息，能否推断出read mapping到的区域对应的参考基因组序列? (提示:参考中对于MD tag的介绍) * （4）软件安装和资源文件的下载也是生物信息学实践中的重要步骤。请自行安装教程中未涉及的[bwa](https://github.com/lh3/bwa)软件，从[UCSC Genome Browser](https://hgdownload.soe.ucsc.edu/downloads.html)下载Yeast (S. cerevisiae, sacCer3)基因组序列。使用`bwa`对Yeast基因组`sacCer3.fa`建立索引，并利用`bwa`将`THA2.fa`，mapping到Yeast参考基因组上，并进一步转化输出得到`THA2-bwa.sam`文件。 {% hint style="info" %} **Help: SAM V.S. BAM** 在 NGS（下一代测序）数据处理中，SAM 和 BAM 文件**本质上存储的是相同的内容**，它们的主要区别在于存储格式和读取效率。 * SAM (Sequence Alignment Map): 是一种纯文本格式。它是人类可读的，你可以直接用 `cat` 或 `head` 命令在终端查看其中的内容。 * BAM (Binary Alignment Map): 是 SAM 文件的二进制压缩版本。它不可直接阅读，主要供计算机程序（如 IGV、Samtools）高效处理。 在实际的生物信息学流程中，我们几乎总是将 SAM 转换为 BAM，原因如下： 1. 节省空间： BAM通常比 SAM **小 75% 左右**，NGS 数据通常有数亿行，SAM 文件动辄几百 GB，转换为 BAM 可以极大地缓解存储压力。 2. 建立索引 (.bai)： 只有 BAM 文件可以建立索引。有了索引，分析软件（如 IGV 基因组浏览器）就可以直接跳转到基因组的特定位置（如第 17 号染色体），而不需要从头开始扫描整个文件。 3. 计算效率： 现代生物信息学工具（如 Samtools, GATK, BWA）都针对二进制流进行了优化。 因此， 推荐mapping生成 bam，例如：`bwa mem sacCer3.fa THA2.fa | samtools sort -o THA2-bwa.sorted.bam` 不推荐生成 sam，例如： `bwa mem sacCer3.fa THA2.fa > THA2-bwa.sam` {% endhint %} ## 7) More Reading * 如上所述，mapping本身用到的都是一些现成的，高度优化的工具，需要我们自己完成的工作不多。 * 很多情况下我们直接拿到的测序数据在mapping前需要进行对各种adaptor和linker的移除、对质量不好的raw reads的trim和丢弃等data clean和QC工作。这是需要一定的经验的，而且对于不同测序方法产生的数据，预处理的方法也有所不同。这需要我们对产生数据的实验方法有清晰了解，illumina团队把现有的所有二代测序建库技术都收集并整理成标准，免费查看。请参考[Sequence Method Explorer](https://www.illumina.com/science/sequencing-method-explorer.html) (中文版链接为)。 * 对于数据的QC和预处理，也有大量工具可以直接使用，最常用的有: * [fastqc](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/): 统计多种QC指标，以html的形式进行可视化 * [cutadapt](https://cutadapt.readthedocs.io/en/stable/): 用于去除adaptor和低质量序列 * [trim\_galore](https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/):对[cutadapt](https://cutadapt.readthedocs.io/en/stable/)进行了封装，自动识别常见adaptor * [fastp](https://github.com/OpenGene/fastp):国人开发的fastq预处理工具 * [fastx\_toolkit](http://hannonlab.cshl.edu/fastx_toolkit/) 一个比较早期的fastq预处理工具 ## 8) Teaching Video * see Videos in the \*\*\*\* [**Files needed**](https://courses.ncrnalab.org/files) ## Take a break **Illumina & Affymetrix** ![Illumina](/files/-LPVvEv4rA7I2GCMcz23) > “Illumina公司的市场数据实在是非常美妙的东西。”拥有个人博客的基因研究人员丹尼尔·麦克阿瑟（Daniel Macarthur）说，“它是如此地纯净，到了令人吃惊的地步。” 当Illumina公司目前的股价净值比高达84倍的时候，高盛仍然建议买入，声称该公司很有可能继续保持其在DNA测序领域里的领导地位。 Illumina毫无疑问是这个时代科技公司的榜样。在生物学仪器制造领域，跟Thermo，Life等拥有几十年历史的公司相比，1998年起源于加州圣地亚哥的Illumina显然还是个年轻小伙。从1999年25个人的规模、130万美元的年营业额，到2013年超过3200人的规模、14.21亿美元的年营业额。 今天，Illumina公司几乎垄断了所有的二代测序（NGS：Next Generation Sequencing）市场。但是，Illumina公司最初是一家生产DNA芯片（microarray chip)的公司，这是一种侦测DNA变异的早期技术。而且，那时的Illumina公司还落后于该市场缔造者Affymetrix公司。 Illumina今天已经赶超了Affymetrix，并将其远远抛之脑后。最主要的原因，在于它从DNA芯片到DNA测序上的成功转型。Illumina公司的成功，很多人归功于其CEO，Jay Flatley，在战略上的敏锐。他很可能称得上是生物科技领域里的最佳CEO。2006年，Flatley说服Illumina公司董事会以6亿美元的价格收购了一家叫Solexa的开发NGS技术的小公司（同类型的NGS技术公司还有454等，因为其技术、成本和市场上的弱势，今天已经很少有人听说了），借此大力押注DNA测序。 如今，illumina在二代测序（NGS：Next Generation Sequencing）乃至整个测序市场占据领导地位。引用illumina官方说法：“世界上90%以上的测序数据都由illumina仪器产生”，不较真的话，这句话确实在某种程度上反应了illumina雄踞NGS市场的现状（下图是illumina 测序产品发布时间线） 。尤其是HiSeq系列测序仪的问世，以通量高，产量大，生产规模著称，能够快速、经济的进行大规模平行测序，在大型全基因组测序，全转录组，全外显子组测序，靶向基因测序方面优势明显。 ![](/files/-LPVvEv69CN0D0_T1u9I) --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://book.ncrnalab.org/teaching/part-iii.-ngs-data-analyses/1.mapping.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.