基因组重测序全流程(简易版)
已知达松维尔拟诺卡氏菌亚种是会导致人类放线菌瘤的环境生物，该样本是来自Keddieii血根杆菌DSM 10542的通用样品，旨通过基因组重测序探索其和参考基因组有何不同，找出基因组变异信息。
1.需要的软件
• 软件名：Aspera 版本号：4.0.2.38
• 软件名：sratoolkit 版本号：2.11.1
• 软件名：FastQC 版本号：0.11.7
• 软件名：Trimmomatic版本号：0.38
• 软件名：bwa 版本号：0.7.17-r1188
• 软件名：samtools 版本号：1.7
• 软件名：Annovar 版本：$Date: 2019-10-24 00:05:27 -0400 (Thu, 24 Oct 2019)
• 软件名：bcftools 版本：1.7
2.数据下载
2.1下载sra数据
1.在NCBI官网的搜索框输入SRR022534，可以看到该菌的目前研究以及进展，在最下面Runs点击SRR022534，在跳转界面点击Data access，会看到一个网址，复制网址；
2.服务器建re_seq文件夹，wget下载亚种基因组文件。
wget https://sra-pub-runodp.s3.amazonaws.com/sra/SRR022534/SRR022534
或者prefetch下载
prefetch SRR022534
也可以选择ascp下载，不过如果基因组文件不是特别大不推荐这种方式，比较麻烦，需要到EBI-ENA数据库赋值ftp地址。
2.2.下载参考基因组序列和基因组注释文件
网址：ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/001/877/055/GCA_001877055.1_ASM187705v1/GCA_001877055.1_ASM187705v1_genomic.fna.gz
1.进入NCBI官网，输入GCA_001877055，在界面右边点击FTP directory for GenBank assembly，在下一个界面复制参考基因组.fna.gz文件和基因组注释文件.gff.gz；
2.wget下载。
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/001/877/055/GCA_001877055.1_ASM187705v1/GCA_001877055.1_ASM187705v1_genomic.fna.gz
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/001/877/055/GCA_001877055.1_ASM187705v1/GCA_001877055.1_ASM187705v1_genomic.gff.gz
3.主要分析步骤和结果
3.1系列与参考基因组的比对
3.1.1数据文件的格式转换
在NCBI数据库检索序列号SRR033534查看experiment得知测试数据类型为454的单末端测序。
fastq-dump SRR022534.sra
3.1.2质量评估
1.fastqc质量评估，文件最好用绝对路径；
touch fastp.sh
vim fastp.sh
#!/bin/bash
for i in 1 2 3 4
do
fastp -i SRR137492${i}.fastq.gz -o SRR137492${i}.fq.gz
done
#Esc : wq! 保存并退出
chmod 777 fastp.sh
./fastp.sh
2.本地查看fastqc report，由结果可知，测序数据质量不是特别高，Per base sequence quality中有很多下四分位小于10或中位数小于25，所以下一步要先数据过滤，去除低质量序列。
3.1.3测序数据过滤
Trimmomatic参数说明
PE：双端测序
-phread33：设置碱基的质量格式
ILLUMINACLIP：
TruSeq3-PE.fa： <fastaWithAdaptersEtc>是fasta格式的接头序列文件路径
2： <seed mismatches>是将接头序列的一段(不超过16bp)作为seed，与reads进行比对，能够容忍的最大错配(mismatch)数，这里是最多2个错配
30： <palindrome clip threshold>是 a-R1, a-R2的比对分值阈值，达到阈值，才进行切除，这里设置阈值为30（大约比对50bp碱基）
10：<simple clip threshold>是任意(接头)序列与read比对最低分阈值，大于这个阈值，才进行切除，这里设置阈值为10（大约比对17bp碱基）
LIDINGWINDOW:5:20：设置滑动窗口阈值，以为5bp为窗口，这5bp碱基的平均质量值低于20，要进行切除
:LEADING:20：设置从reads起始开始，去除质量低于阈值或为’N’的碱基，直到达到阈值不再去除，这里设置阈值为20
TRAILING:20：设置从reads末尾开始，去除质量低于阈值的碱基或为’N’的碱基直到达到阈值不再去除，这里设置阈值为3，这种方法是去除特定的illumina平台低质量区域(由于illumina会将低质量碱基标记为2)
MINLEN:75 设置read切除后的最短长度，这里设置长度至少为36bp，长度小于36bp的reads被去除
fastqc /disk1/202031107010114/re_seq/trim_out/SRR022534.fastq.gz
3.1.4fastqc重新质量评估
fastqc /disk1/202031107010114/re_seq/trim_out/SRR022534.fastq.gz
3.1.5建立参考基因组索引
gunzip disk1/202031107010114/re_seq/GCA_001877055.1_ASM187705v1_genomic.fna.gz
bwa index disk1/202031107010114/re_seq/GCA_001877055.1_ASM187705v1_genomic.fna
3.1.6测序数据比对到参考基因组得到sam文件
bwa mem disk1/202031107010114/re_seq/GCA_001877055.1_ASM187705v1_genomic.fna disk1/202031107010114/re_seq/trim_out/SRR022534_out.fastq.gz >bwa_mem_SRR022534.sam
3.1.7sam文件转bam文件
samtools faidx disk1/202031107010114/re_seq/GCA_001877055.1_ASM187705v1_genomic.fna
samtools view -bhS -t GCA_001877055.1_ASM187705v1_genomic.fna.fai -o bwa_mem_SRR022534.bam bwa_mem_SRR022534.sam
3.1.8为bam文件排序
samtools sort bwa_mem_SRR022534.bam -o bwa_mem_SRR022534.sorted.bam
3.1.9为bam文件建立索引 显示基因组比对情况
samtools index bwa_mem_SRR022534.sorted.bam
samtools index bwa_mem_SRR022534.sorted.bam
3.1.10测试每个基因组每个位点的测试深度并统计比对结果
samtools depth bwa_mem_SRR022534.sorted.bam >>depth.txt
samtools flagstat bwa_mem_SRR022534.sorted.bam
3.2变异点检测
3.2.1去除pcr重复
samtools rmdup bwa_mem_SRR022534.sorted.bam bwa_mem_SRR022534_nopcr.bam
3.2.2生成bcf文件
samtools mpileup -gf GCA_001877055.1_ASM187705v1_genomic.fna bwa_mem_SRR022534_nopcr.bam >bwa_mem_SRR022534.bcf
3.2.3基因变异检测和筛选
bcftools call -vm bwa_mem_SRR022534.bcf -o bwa_mem_SRR022534.variants.bcf
bcftools view -v snps,indels bwa_mem_SRR022534.variants.bcf > bwa_mem_SRR022534.snps.vcf
3.3变异位点注释
3.3.1生成annovar输入文件
convert2annovar.pl -format vcf4 bwa_mem_SRR022534.snps.vcf > bwa_mem_SRR022534.snps.avinput
3.3.2自定义注释数据框
gunzip disk1/202031107010114/re_seq/GCA_001877055.1_ASM187705v1_genomic.gff.gz
wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/gff3ToGenePred
chmod 777 gff3ToGenePred
./gff3ToGenePred -useName GCA_001877055.1_ASM187705v1_genomic.gff 7055-genome_refGene.txt
cut -f 12 7055-genome_refGene.txt >column1.txt
cut -f 2-15 7055-genome_refGene.txt >column_else.txt
paste column1.txt column_else.txt >7055-genome_new_refGene.txt
retrieve_seq_from_fasta.pl -format refGene -seqfile GCA_001877055.1_ASM187705v1_genomic.fna -outfile 7055-genome_new_refGeneMrna.fa 7055-genome_
cp 7055-genome_new_refGene* ~/Biosofts/annovar/humandb/
3.3.3注释变异基因位点
annotate_variation.pl --geneanno --dbtype refGene --buildver 7055-genome_new bwa_mem_SRR022534.snps.avinput ~/Biosofts/annovar/humandb/
4.批量处理的脚本
因为我所分析的sra文件只有一个，无法更好地说明批量处理过程，就找了一组RNA-Seq数据，二者前面的分析流程基本一致。数据来自研究：抗体富集前后培养物中伯氏疏螺旋体基因表达中富含aBa的体外培养Bb代表SRR22149531和体外培养的Bb代表SRR22149536。
4.1prefetch批量下载sra文件
touch prefetch.sh
vim prefetch.sh
#!/bin/bash
for i in 1 6
do
prefetch SRR2214953${i}
done
#Esc : wq! 保存并退出
chmod 777 prefetch.sh
./prefetch.sh
4.2fastq-dump批量从sra文件提取fastq文件
touch fastq_dump.sh
vim fastq_dump.sh
#!/bin/bash
for i in SRR*
do
echo ${i}#打印输出文件名
fastq-dump --gzip ${i}/${i}.sra
#可根据数据需要选择--spilt-files/--spilt-3/--spilt-spot
done
#Esc : wq! 保存并退出
chmod 777 fastq_dump.sh
./fastq_dump.sh
4.3fastqc批量质控
touch fastqc.sh
vim fastqc.sh
#！/bin/bash
for i in 1 6
do
fastqc SRR2214953${i}
done
#Esc : wq! 保存并退出
chmod 777 fastqc.sh
4.4bwa批量比对并生成排好序的bam文件
touch bwa_samtools.sh
vim bwa_samtools.sh
#!/bin/bash
for i in 1 6;do
{
bwa index GCA_000181575.2_ASM18157v2_genomic.fna;
bwa mem GCA_000181575.2_ASM18157v2_genomic.fna SRR2214953${i}.fastq.gz >bwa_mem_SRR2214953${i}.sam;
samtools faidx GCA_000181575.2_ASM18157v2_genomic.fna;
samtools view -bhS -t GCA_000181575.2_ASM18157v2_genomic.fna.fai -o bwa_mem_SRR2214953${i}.bam bwa_mem_SRR2214953${i}.sam;
samtools sort bwa_mem_SRR2214953${i}.bam -o bwa_mem_SRR2214953${i}.sorted.bam;}
done
#Esc : wq! 保存并退出
chmod 777 bwa_samtools.sh
./bwa_samtools.sh
基因组重测序全流程(简易版)相关推荐
- 一个全基因组重测序分析实战
Original 2017-06-08 曾健明 生信技能树 这里选取的是 GATK best practice 是目前认可度最高的全基因组重测序分析流程,尤其适用于 人类研究. PS:其实本文应该属于 ...
- signature=d363d26bda212f777fef81d270ecd42b,基于DNA-pooling全基因组重测序初步筛查CAD易感基因变异位点...
摘要: 目的:应用DNApooling全基因组重测序技术初步筛查CAD易感基因变异位点并进行功能分析.方法:分别收集CAD病例组血液样本和正常对照组血液样本各100例.提取病例组和对照组DNA,制作D ...
- 全基因组重测序基础及高级分析知识汇总
全基因组重测序基础及高级分析知识汇总 oddxix 已关注 2018.09.20 17:04 字数 11355 阅读 212评论 0喜欢 6 转自:http://www.360doc.com/cont ...
- Flutter IOS 新建打包发布全流程 2023 版
大家好,我是 17. 上一篇写完 Flutter Android 打包保姆式全流程 2023 版 后,小伙伴 MannaYang 和 Mapleeeeeee 留言说要看 IOS 打包的流程.于是 17 ...
- Flutter Android 打包保姆式全流程 2023 版
大家好,我是 17. Flutter 打包的文章一共有两篇 Flutter Android 打包保姆式全流程 2023 版 Flutter IOS 新建打包发布全流程 2023 版 本篇介绍 Andr ...
- 全基因组重测序揭示了野生大豆的局部适应和分化的特征
近日,凌恩生物客户南京农业大学在<Evolutionary Applications>发表题为"Whole-genome resequencing reveals signatu ...
- 中国人泛基因组发了 Nature，基因组重测序的CNS值得拥有
33333福利公告:前 10期<临床基因组/外显组数据分析实战>线上/线下课程已圆满结束.现于2023年6月30~7月1日,在北京安排第十一期课程. 线上课是通过腾讯会议实时直播线下课,实 ...
- 全基因组重测序基础分析
本文仅为个人记录 流程说明 从零开始完整学习全基因组测序(WGS)数据分析:第3节 数据质控 2015a 2016a 全基因组测序数据获取后应该怎么分析 samtools samtools常用命令详解 ...
- 2020.8.5丨细菌基因组二代测序组装流程梳理
目录 step:1 准备工作 step:2 数据过滤 比对参考基因组 安装BWA 建立索引 bwa一共有三种比对算法,这里使用mem 数据过滤 安装samtools 筛选三步走 step:3 组装 s ...
最新文章
- 基于Flink的在线机器学习系统架构探讨
- Windows和Linux系统下，虚拟环境安装的全面说明和详细步骤
- 通过PRINT过程制作报表
- Android之ViewPager讲解
- 四种排序（冒泡、插入、递归、选择）
- Selenium3自动化测试——11. 下拉框处理
- 对一道面试题的总结与扩展思考（关于一笔画问题的数学分析）
- 阿里云有一群 “猪猪侠”
- aws集群重启_使用自动伸缩组在AWS中运行安全数据库集群
- leetcode 303. 区域和检索 - 数组不可变
- Hibernate 多表关联
- iOS 13.2正式版放出 支持AirPodsPro
- C++_选择结构_循环结构_for循环_敲桌子案例_嵌套循环_乘法口诀案例_跳转语句break---C++语言工作笔记018
- 虚拟主机如何创建svn服务器,虚拟主机搭建svn
- Linux中 查看mysql配置文件位置
- android中pdf转换成图片格式,Android-PDF转图片
- elasticsearch 版本区别
- 1.1【气宇轩昂】《踏雪》
- 羽绒枕头行业调研报告 - 市场现状分析与发展前景预测（2021-2027年）
- 关于H.265/HEVC视频压缩标准相较H.264/AVC节省50%左右的带宽方案的推荐