全基因组测序数据分析---WGS主流程

最新推荐文章于 2024-04-23 20:56:28 发布
最新推荐文章于 2024-04-23 20:56:28 发布
阅读量4.3k 收藏 21
点赞数 1
文章标签： 数据分析 大数据
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/dengshaoxiong/article/details/127430403
版权

WGS数据分析流程

1 准备阶段

部署好相关的软件和工具
BWA (Burrow-Wheeler Aligner) Version 0.7.17-r1188 解压、编译
Samtools Version: 1.16.1

解压 tar jxvf samtools-1.16.1.tar.gz
进入目录 cd samtools-1.16.1
配置 ./configure --prefix=~/biosoft/samtools-1.16.1
编译安装 make
make install

Picard 直接下载java包picard.jar
GATK gatk-4.3.0.0 下载后不用编译直接使用
上述软件均需添加到环境变量

2 数据预处理

2.1 构建索引

参考基因组索引的构建

bwa index Sus_scrofa.Sscrofa11.1.dna.toplevel.fa

java -jar /home/dengsx/biosoft/picard.jar CreateSequenceDictionary \
R=/home/dengsx/publicdata/reference_genome/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa \
O=/home/dengsx/publicdata/reference_genome/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.dict

samtools faidx Sus_scrofa.Sscrofa11.1.dna.toplevel.fa ##fai索引的构建

后续GATK的使用需要用到多种索引类型，要用多个软件创建

dict索引文件也可以通过gatk来获取

gatk CreateSequenceDictionary -R E.coli_K12_MG1655.fa -O E.coli_K12_MG1655.dict && echo "** dict done **"(本人没试过)
###需要注意的是.dict文件的名字前缀需要和fasta的一样，并跟它在同一个路径下，这样GATK才能够找到。

dbSNP索引的构建（2.6 BQSR之前做好就行）

java -jar /home/dengsx/biosoft/gatk-4.3.0.0/gatk-package-4.3.0.0-local.jar IndexFeatureFile --input /home/dengsx/publicdata/dbsnp/sus_scrofa.nospace.vcf > /home/dengsx/publicdata/dbsnp/sus_scrofa.nospace.vcf.index.log 2>&1

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-hErFUiEQ-1666255841603)(Snipaste_2022-10-19_21-52-02.png)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-nuH0BtQJ-1666255841604)(Snipaste_2022-10-19_21-53-48.png)]
这一步完成之后，我们就可以将read比对至参考基因组了

2.2 bwa比对

将比对的输出结果直接重定向到一份*.sam文件中，这类文件是BWA比对的标准输出文件，。但SAM文件是文本文件，一般整个文件都非常巨大，因此，为了有效节省磁盘空间，用samtools将它转化为BAM文件（SAM的特殊二进制格式），而且BAM会更加方便于后续的分析。

bwa mem -t 5 -R '@RG\tID:YF62_E7\tPL:UNKNOWN\tLB:library2\tSM:YF62_E7' /home/dengsx/publicdata/reference_genome/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa YF62_E7_1_clean.fq.gz YF62_E7_2_clean.fq.gz | samtools view -S -b - > YF62_E7.bam

-t，线程数，我们在这里使用4个线程；-R 接的是Read Group的字符串信息，这是一个非常重要的信息，以@RG开头，它是用来将比对的read进行分组的。不同的组之间测序过程被认为是相互独立的，这个信息对于我们后续对比对数据进行错误率分析和Mark duplicate时非常重要。

在Read Group中，有如下几个信息非常重要：

  • ID，这是Read Group的分组ID，一般设置为测序的lane ID（不同lane之间的测序过程认为是独立的），下机数据中我们都能看到这个信息的，一般都是包含在fastq的文件名中；
  • SM，sample_name 同一个个体的sampleID必须一样，有时候我们测序的数据比较多的时候，那么可能会分成多个不同的lane分布测出来，这个时候SM名字就是可以用于区分这些样本；
  • PL，指的是所用的测序平台，这个信息不要随便写！在GATK中，PL只允许被设置为：ILLUMINA,SLX,SOLEXA,SOLID,454,LS454,COMPLETE,PACBIO,IONTORRENT,CAPILLARY,HELICOS或UNKNOWN这几个信息，如果实在不知道，那么必须设置为UNKNOWN；
  • LB，测序文库的名字，这个重要性稍微低一些，主要也是为了协助区分不同的group而存在。文库名字一般可以在下机的fq文件名中找到，如果上面的lane ID足够用于区分的话，也可以不用设置LB。
    以上，我们就完成了read比对的步骤。接下来是排序：

2.3 merge个体

最好是在这一步就merge，后面的操作都是针对merge后的文件

samtools merge -o srr.bam srr1.bam srr2.bam

2.4 排序

samtools sort -@ 1 -m 64G -O bam -o A65.sorted.bam A65.bam

其中，-@，用于设定排序时的线程数；-m，限制排序时最大的内存消耗，-O 指定输出为bam格式；-o 是输出文件的名字。建议在做类似分析的时候在文件名字将所做的关键操作包含进去，因为这样即使过了很长时间，当你再去看这个文件的时候也能够立刻知道当时对它做了什么。

2.5 去除重复序列（或者标记重复序列）

数据预处理

for i in `ls *.sorted.bam`
do
java -jar /home/dengsx/biosoft/picard.jar MarkDuplicates \
REMOVE_DUPLICATES=false \
I=${input}/${i} \
O=${output}/${i%\.bam}.markup.bam \
M=${output}/${i%\.bam}.markup_metrics.txt
  1. 这里只把重复序列在输出的新结果中标记出来，但不删除。如果我们非要把这些序列完全删除的话可以设置参数REMOVE_DUPLICATES=true建议使用第一种做法，只是标记出来，并留存这些序列，以便在你需要的时候还可以对其做分析。
  2. 这一步完成之后，我们需要为sample_name.sorted.markdup.bam创建索引文件，它的作用能够让我们可以随机访问这个文件中的任意位置，而且后面的“局部重比对”步骤也要求这个BAM文件一定要有索引，命令如下：
cd /home/dengsx/publicdata/markup_out
for markup_bam in `ls *.sorted.markup.bam`
do
samtools index ${markup_bam}
done

在重新校正碱基质量值（BQSR）之前把相同个体的bam文件merge，将同个样本的所有比对结果合并成唯一一个大的BAM文件

samtools merge -o srr.sorted.markup.bam srr1.sorted.markup.bam srr2.sorted.markup.bam

2.6 重新校正碱基质量值（BQSR）

主要是通过机器学习的方法构建测序碱基的错误率模型，然后对这些碱基的质量值进行相应的调整。
这里包含了两个步骤：

  • 第一步，BaseRecalibrator，这里计算出了所有需要进行重校正的read和特征值，然后把这些信息输出为一份校准表文件（sample_name.recal_data.table）

  • 第二步，ApplyBQSR ，这一步利用第一步得到的校准表文件（sample_name.recal_data.table）重新调整原来BAM文件中的碱基质量值，并使用这个新的质量值重新输出一份新的BAM文件。

数据准备
dbSNP
数据处理

  1. 解压sus_scrofa.vcf.gz文件
  2. the dbsnp vcf file contains axiom snp array snps that offends the gatk
    because they have white space in the INFO field, remove these (去除空白格)
awk -F "\t" '!($8 ~ /\s/)' /home/dengsx/publicdata/dbsnp/sus_scrofa.vcf > /home/dengsx/publicdata/dbsnp/sus_scrofa.nospace.vcf
  1. 构建.idx索引
java -jar /home/dengsx/biosoft/gatk-4.3.0.0/gatk-package-4.3.0.0-local.jar IndexFeatureFile --input /home/dengsx/publicdata/dbsnp/sus_scrofa.nospace.vcf > /home/dengsx/publicdata/dbsnp/sus_scrofa.nospace.vcf.index.log 2>&1

BQRS重新校正碱基质量值

BQRS 第一步(BaseRecalibrator)

fasta=/home/dengsx/publicdata/reference_genome/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa
input=/home/dengsx/publicdata/markup_out/

gatk --java-options "-Xmx10G -Djava.io.tmpdir=./" BaseRecalibrator \
    -R ${fasta} \
    -I ${input}/A65.sorted.markup.bam \
    --known-sites /home/dengsx/publicdata/dbsnp/sus_scrofa.nospace.vcf \
    -O A65.recal_data.table
  • 这里计算出了所有需要进行重校正的read和特征值，然后把这些信息输出为一份校准表文件（sample_name.recal_data.table）

# 不能只用bwa构建的索引来做BQSR
# 需要的索引类型很多(参考上面索引的构建)

  • 参考基因组的索引.fai、.dict
java -jar /home/dengsx/biosoft/picard.jar CreateSequenceDictionary R=/home/dengsx/publicdata/reference_genome/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa O=/home/dengsx/publicdata/reference_genome/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.dict
samtools faidx Sus_scrofa.Sscrofa11.1.dna.toplevel.fa
  • sus_scrofa.nospace.vcf的索引.idx
  • bam文件的索引

BQRS 第二步(ApplyBQSR)

fasta=/home/dengsx/publicdata/reference_genome/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa
input=/home/dengsx/publicdata/markup_out/
bqsr=/home/dengsx/publicdata/BQSR/
gatk --java-options "-Xmx10G -Djava.io.tmpdir=./" ApplyBQSR \
    -R ${fasta} \
    -I ${input}/A65.sorted.markup.bam \
    --bqsr-recal-file ${bqsr}/A65.recal_data.table \
    -O A65.sorted.markdup.BQSR.bam
  • 这一步利用第一步得到的校准表文件（sample_name.recal_data.table）重新调整原来BAM文件中的碱基质量值，并使用这个新的质量值重新输出一份新的BAM文件。

3 变异检测

变异检测流程图
我们这里使用GATK HaplotypeCaller模块对样本中的变异进行检测，它也是目前最适合用于对二倍体基因组进行变异（SNP+Indel）检测的算法。

一般来说，在实际的WGS流程中对HaplotypeCaller的应用有两种做法，差别只在于要不要在中间生成一个gVCF：

  • 第一种，直接进行HaplotypeCaller，这适合于单样本，或者那种固定样本数量的情况，也就是执行一次HaplotypeCaller之后就老死不相往来了。否则你会碰到仅仅只是增加一个样本就得重新运行这个HaplotypeCaller的坑爹情况（即，N+1难题），而这个时候算法需要重新去读取所有人的BAM文件，这将会是一个很费时间的痛苦过程；
  • 第二种，每个样本先各自生成gVCF，然后再进行群体joint-genotype。这其实就是GATK团队为了解决（1）中的N+1难题而设计出来的模式。gVCF全称是genome VCF，是每个样本用于变异检测的中间文件，格式类似于VCF，它把joint-genotype过程中所需的所有信息都记录在这里面，文件无论是大小还是数据量都远远小于原来的BAM文件。这样一旦新增加样本也不需要再重新去读取所有人的BAM文件了，只需为新样本生成一份gVCF，然后重新执行这个joint-genotype就行了。

基因组上各个不同的染色体之间其实是可以理解为相互独立的（结构性变异除外），也就是说，为了提高效率我们可以按照染色体一条条来独立执行这个步骤，最后再把结果合并起来就好了，这样的话就能够节省很多的时间

input=/home/dengsx/publicdata/BQSR
output=/home/dengsx/publicdata/haplotypeCaller
reference=/home/dengsx/publicdata/reference_genome/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa
cd /home/dengsx/publicdata/haplotypeCaller
ls ../BQSR/*.sorted.markdup.BQSR.bam | while read line
do
gatk --java-options "-Xmx10G -Djava.io.tmpdir=./" HaplotypeCaller  \
    -I $input/${line##../BQSR/} \
    -R $reference \
    -ERC GVCF \
    -ploidy 2 \
    -O $output/${line##../BQSR/}.g.vcf > /home/dengsx/publicdata/haplotypeCaller/haplotypeCaller.log 2>&1
done
wait # this is necessary because both processes need to complete for the outside call to check on logs
echo $(date) done.main.process

参考：碱基矿工

dengshaoxiong
关注
  • 1
    点赞
  • 踩
  • 21
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
wgs-Analysis-process:使用gatk进行wgs全基因组分析寻找SNP变异的流程
04-29
wgs全基因组序列比对流程 用到的软件 过程步骤 一. 下载准备需要的文件 下载参考序列基因组文件 1.建立索引 bwa index ref.fasta 完成之后 会看到几个ref.fasta为前缀的文件 为参考序列生成dict文件 gatk CreateSequenceDictionary -R ref.fasta -O ref.dict samtools 建索引 samtools faidx ref.fasta 下载测序文件 fastaq-dump --split-files SRR***** 下载的文件是双末端测序从两端读的read1和read2 >> 用bgzip压缩 bgzip seq1_.fasta bgzip seq2_.fasta 二.处理文件 将read比对到参考基因组 bwa mem -t 4 -R '@RG\tID:foo\tPL:illumina\tSM:
WGS-全基因组测序
weixin_57590035的博客
03-01 286
过程：提取基因组DNA→构建测序文库（DNA片段化、两端加接头）→高通量测序→测序数据预处理（质控）→参考序列比对（read分布统计、测序覆盖度统计、测序深度统计）→分析（基因变异分析、SNP分析、InDel分析、群体分析、疾病相关分析）流程：原始数据质控→过滤低质量reads→read比对→排序sort→去重复→布局重比对→碱基质量重矫正BQSR→变异检测→(Mergeloptional→joint Genotype)→变异质控和过滤VQSR。优势：snp、插入缺失突变。
《全基因组测序WGS数据分析——1.DNA测序技术》
wx1871428的博客
07-05 2205
WGS(Whole Genome Sequencing) 指将物种细胞里面完整的基因组序列全部DNA，检测并排列，此技术几乎能够鉴定出基因组上任何类型的突变。 对于人类来说，全基因组测序的价值是极大的，它的信息包含了所有基因和生命特征之间的内在关联性，当然也意味着更大的数据解读和更高的技术挑战。 测序，简单来说就是将DNA化学信号转变为计算机可处理的数字信号。 第一代测序技术——sanger法 Sanger法是基于DNA合成反应的测序技术，又称为SBS法、末端终止法。1975年由Sange
WGS（重测序）分析详解与脚本
小猪打呼噜
05-13 1万+
这篇文章很长，超过1万字，是本系列中最重要的一篇，因为我并非只是在简单地告诉大家几条硬邦邦的操作命令。对于新手而言不建议碎片时间阅读，对于有一定经验的老手来说，相信依然可以有所收获。在开始之前，我想先说一句：流程的具体形式其实是次要的，WGS本质上只是一个技术手段，重要的是，我们要明白自己所要解决的问题是什么，所希望获取的结果是什么，然后再选择合适的技术。这是许多人经常忽视的一个重要问题。 好了，以下进入正文。 这是WGS数据分析的流程图。流程的目的是准确检测出每个样本（这里特指人）基因组中的变异.
GWAs——全基因组关联分析流程
SUN5_The_answer的博客
09-19 7200
GWAs基本流程
【生信】全基因组测序（WGS）
热门推荐
weixin_40695088的博客
03-07 3万+
1、全基因组测序(WGS) 的定义 2、GWS流程 2.1准备工作——分析软件 2.2原始数据质控 2.3数据预处理 2.4变异检测
WGS完整流程介绍(原始数据质控、数据预处理、变异检测、数据注释)
qq_36717487的博客
03-10 1万+
111
WGS数据分析流程学习与开发过程全纪录（1）
littleblack201608的博客
02-14 4728
A. 开发环境及数据准备工作 A1. 开发环境搭建 安装miniconda3，然后在其中创建WGS的开发环境： # 1. Install miniconda3. wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh . bash Miniconda3-latest...
1.全基因组重测序数据分析介绍
weixin_53682198的博客
04-14 878
二代测序原理及fastq数据解读、全基因组重测序数据分析内容
WGS数据分析流程学习与开发过程全纪录（2）
littleblack201608的博客
02-14 1361
上一篇主要是跟着别的教程走了一遍大肠杆菌的数据分析流程，这一篇开始就一块一块地具体学习。 本篇主要学习的是数据质控。 1. FastQC a. 下载并安装、配置FastQC (wgs) [xuquan@zhswlcserver:conda]$cd ~/conda/softwares/ (wgs) [xuquan@zhswlcserver:softwares]$wget http://www.bio...
基因组二代测序数据的自动化分析流程
05-10
基因组二代测序数据的自动化分析流程基因组二代测序数据的自动化分析流程
人工智能-数据分析-川牛膝全基因组高通量测序及初步数据分析.pdf
07-01
人工智能-数据分析-川牛膝全基因组高通量测序及初步数据分析.pdf
pyrus-seq:Pyrus 分析双端全基因组测序数据-开源
06-28
Pyrus 分析双端和非双端全基因组测序数据。 该程序旨在使用源自数据集本身的参数来识别和统计数据中识别出的潜在染色体重排。
A_rabiei_WGS_analysis:狂犬病菌分离株的全基因组测序分析流程
05-15
2017年，从21种狂犬病壁球菌菌株中提取DNA，并在Illumina HiSeq2500上进行全基因组测序（WGS），产生100 bp短配对末端读数（Macrogen，Korea）。 在第二年（2018年），又从20种其他狂犬病菌株中提取了DNA，并将其...
人工智能-数据分析-基于二代测序的转录组数据分析方法的比较研究.pdf
07-02
人工智能-数据分析-基于二代测序的转录组数据分析方法的比较研究.pdf
数据分析场景，连号相关业务
weixin_50235904的博客
04-20 323
演示的为db2数据库，需要含有窗口函数，或者可以获取到当前数据偏移的上一位数据。
DaPy：实现数据分析与处理
最新发布
爱编程的鱼的博客
04-23 896
DaPy是一个基于Python的开源库，专注于数据分析和处理。它提供了一套简洁而灵活的工具和函数，使开发者能够对数据进行各种操作，如数据清洗、转换、筛选和聚合等。DaPy库的设计目标是帮助开发者在数据分析过程中高效地处理和分析数据，从而得出准确的结论和洞察。DaPy是一个功能强大的Python库，专注于数据分析和处理。它提供了丰富的工具和函数，使开发者能够高效地进行数据清洗、转换、筛选和聚合等操作。通过使用DaPy库，开发者可以提升数据分析的效率和准确性，从而得出准确的结论和洞察。
人工智能对数据分析领域的影响是什么？
WhiteCattle_DATA的博客
04-23 392
在数据分析过程中，数据清洗是一个至关重要的环节，旨在消除数据中的错误、异常和冗余，确保数据质量。而人工智能算法，特别是机器学习和深度学习模型，能够自动地从海量数据中识别模式、预测趋势并做出决策，大大提高了分析的速度和准确性。这一领域的技术革新和应用拓展正在推动数据分析从传统的基于规则的方法转变为基于机器学习和深度学习的自动化方法。总的来说人工智能对数据分析领域的影响是全方位的，从提高分析效率和准确性到促进数据清洗和可视化，再到解决挑战和限制，都在推动着数据分析领域的不断进步和发展。
数据分析-pandas1
别追我我有止咳糖浆
04-23 632
Pandas 是 Wes McKinney 在2008年开发的一个强大的的工具集。Pandas 以 NumPy 为基础（实现数据存储和运算），提供了专门用于数据分析的类型、方法和函数，对数据分析和数据挖掘提供了很好的支持；同时 pandas 还可以跟数据可视化工具 matplotlib 很好的整合在一起，非常轻松愉快的实现数据可视化呈现。Pandas 核心的数据类型是Series（数据系列）、DataFrame（数据窗/数据框），分别用于处理一维和二维的数据，除此之外，还有一个名为Index。
二代测序数据分析流程
07-28
二代测序数据分析流程包括以下几个步骤： 1. 数据质控：对测序数据进行质量评估和过滤，去除低质量的reads和污染序列，以确保后续分析的准确性和可靠性。 2. 序列比对：将过滤后的reads与参考基因组或转录组进行比对，以确定每个read的起始位置和方向。 3. 变异检测：通过比对结果，检测样本中的单核苷酸变异（SNVs）、插入缺失（indels）等遗传变异。 4. 基因表达分析：根据比对结果，计算每个基因的表达水平，包括基因的读数、FPKM（每百万读数的碱基数）或TPM（每百万转录本的碱基数）等。 5. 基因差异表达分析：比较不同条件下的基因表达水平，识别差异表达的基因，并进行统计学分析和功能富集分析。 6. 功能注释：对检测到的变异和差异表达基因进行功能注释，包括基因本体（Gene Ontology）分析、通路富集分析等，以了解其生物学功能和相关通路。 7. 数据可视化：将分析结果以图表、热图、散点图等形式进行可视化展示，以便更好地理解和解释数据。 需要注意的是，二代测序数据分析流程可能因具体研究目的和数据类型的不同而有所差异，上述步骤仅为一般流程的概述。具体的数据分析流程还需要根据实际情况进行调整和优化。\[1\]\[2\]\[3\] #### 引用[.reference_title] - *1* *2* *3* [illumina 二代测序原理及过程](https://blog.csdn.net/zea408497299/article/details/124957981)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v91^insert_down1,239^v3^insert_chatgpt"}} ] [.reference_item] [ .reference_list ]

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

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

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值