liuyujie0136's Website

Logo

A website for self learning, collecting and sharing.

Contact Me

生物信息学常用Linux平台软件简介

Aspera

NCBI-SRA和EBI-ENA数据库

SRA数据库(Sequence Read Archive)

隶属NCBI (National Center for Biotechnology Information)，它是一个保存大规模平行测序原始数据以及比对信息和元数据 (metadata) 的数据库，所有已发表的文献中高通量测序数据基本都上传至此，方便其他研究者下载及再研究。其中的数据则是通过压缩后以.sra文件格式来保存的，SRA数据库可以用于搜索和展示SRA项目数据，包括SRA主页和 Entrez system，由 NCBI 负责维护。SRA数据库中的数据分为Studies, Experiments, Samples和相应的Runs四个层次：

ENA数据库(European Nucleotide Archive)

隶属EBI (European Bioinformatics Institute)，功能等同SRA，并且对保存的数据做了注释，界面相对于SRA更友好，对于有数据需求的研究人员来说，ENA数据库最诱人的点应该是可以直接下载fastq.gz文件，且可以方便地确认下载数据的完整性。

SRA/FASTQ文件下载方式

需要获取他人发表的公开测序数据，来帮助自己的研究领域，下载.sra文件是为了获取该sra相对应的fastq或者sam文件，通过文件格式转换就可以和自己的pipeline对接上，用于直接分析，所以：

  1. 我们需要到SRA或者ENA上搜索我们选择好的SRR号或者SRS号或者SRP号，先在ENA上搜索，如没有再去SRA上搜索，因为ENA下载比SRA快。
  2. 下载数据，从SRA数据库下载数据有多种方法。可以用ascp快速的来下载sra文件，也可以用wget或curl等传统命令从FTP服务器上下载sra文件（但是wget和curl下载的sra文件有时候会不完整），另外NCBI的sratoolkit工具集中的prefetch、fastq-dump和sam-dump也支持直接下载。

使用Aspera从EBI-ENA下载FASTQ数据

推荐使用conda进行安装，使用以下命令：

conda install aspera-cli
# 或：
conda install -c hcc aspera-cli

安装完成后，先了解ascp命令的使用方法与常用参数：

ascp [参数] 目标文件 保存路径

OPTIONS：
-v  verbose mode 唠叨模式，能让你实时知道程序在干啥，方便查错。
-T  取消加密，否则有时候数据下载不了
-i  提供私钥文件的地址，免密从ENA下载，不能少，如：~/miniconda3/pkgs/aspera-cli-3.7.7-0/etc/asperaweb_id_dsa.openssh
-l  设置最大传输速度，一般200m到500m，如果不设置反而速度会比较低
-k  重要！！设置断点续传，一般设为1
-Q  用于自适应流量控制，磁盘限制所需
-P  用于SSH身份验证的TCP端口，一般是33001

下载示例：

批量下载：

bgzip-tabix

tabix可以对NGS分析中常见格式的文件建立索引，从而加快访问速度，不仅支持VCF文件，还支持BED, GFF，SAM等格式。

对于VCF文件，由于SNP位点多，其文件也非常大，需要进行压缩。bgzip可以压缩VCF文件，用法如下：

bgzip example.vcf

压缩之后，原本的example.vcf文件就变成了example.vcf.gz文件。压缩后缀为.gz, 若要解压缩，可以：

bgzip -d example.vcf.gz
gunzip example.vcf.gz

注意：在对VCF文件压缩时，不可以使用gzip来代替bgzip。

对于大型的VCF文件而言，如何快速访问其中的记录也是个难点。tabix可以对VCF文件构建索引，如下：

tabix -p vcf example.vcf.gz

注意输入的VCF文件必须是使用bgzip压缩之后的VCF文件，生成的索引文件为example.vcf.gz.tbi。

构建好索引之后，可以快速的获取指定区域的记录，示例如下：

# 获取位于11号染色体的SNP位点
tabix example.vcf.gz 11
# 获取位于11号染色体上突变位置大于或者等于2343545的SNP位点
tabix example.vcf.gz 11:2343545
# 获取位于11号染色体上突变位置介于2343540到2343596的SNP位点
tabix example.vcf.gz 11:2343540-2343596

snpEff-SnpSift

SnpEff Manual

在进行基因组学研究中，需要对一些SNP和Indel进行分析，而注释它们就需要用到SnpEff，SnpEff是一个变异注释和预测变异影响(氨基酸和核酸改变)的工具。输入文件通常是VCF(variant call format)文件，也接受BED文件(如来自ChIP-Seq peaker)。运行后会输出变异的注释并预测对已知基因的影响。详见官网帮助文档。

示例1: 先用SnpSift选择非低质量位点，再用snpEff进行注释

#!/bin/bash

SnpSift filter \
    "(na FILTER) | (FILTER = 'PASS')" \
    Sample1.snp.vcf.gz | \
    snpEff ann -v \
        -s Sample1.snpEff.summary.html \
        Oryza_sativa \
        > Sample1.snpEff.vcf

示例2: 注释完成后再用SnpSift进一步筛选

#!/bin/bash

SnpSift filter \
    "(ANN[*].IMPACT = 'HIGH') | (ANN[*].IMPACT = 'MODERATE')" \
    Sample1.snpEff.vcf \
    > Sample1.snpEff.HIGH-MODERATE.vcf

# IMPACT = {HIGH, MODERATE, LOW, MODIFIER}
# Other example: "ANN[*].EFFECT has 'missense_variant'"

SnpSift extractFields \
    Sample1.snpEff.HIGH-MODERATE.vcf \
    CHROM POS REF ALT "ANN[0].EFFECT" "ANN[0].GENEID"\
    > Sample1.snpEff.HIGH-MODERATE.gene.txt

cat Sample1.snpEff.summary.genes.txt | sed 's/\t/,/g' > Sample1.snpEff.summary.genes.csv

bedtools

BEDtools Tutorial

BEDTools是可用于genomic features的比较，相关操作及进行注释的工具。而genomic features通常使用Browser Extensible Data (BED) 或者 General Feature Format (GFF)文件表示，用UCSC Genome Browser进行可视化比较。bedtools总共有二三十个工具/命令来处理基因组数据。 比较典型而且常用的功能举例如下：

bedtools genomecov

genomecov可以把alignment的结果文件转为bedgraph格式文件，参考。

有两种方法可以调用之：

bedtools genomecov [OPTIONS] [-i|-ibam] -g (iff. -i)
genomeCoverageBed  [OPTIONS] [-i|-ibam] -g (iff. -i)

这个命令本身并不是设计来做格式转换的，bam2bedgraph也只是其中的一个小功能而已，需要加上-bg参数，就可以report depth in BedGraph format.

示例：

bedtools genomecov -bg -i Sample1.tagAlign -g genome.len > Sample1.bedGraph
bedtools genomecov -bg -ibam Sample2.unique.sorted.bam > Sample2.bedGraph

注意：alignment的文件必须是sorted的，如果是bed格式的比对文件，用-i参数来指定输入文件，需要加入参考基因组的染色体大小记录文件(genome.len)，如果是bam格式的比对文件，用-ibam指定输入文件，而且不需要参考基因组的染色体大小记录文件。

bedtools multicov

multicov可以对RNA-seq比对到各个基因的reads进行计数。

示例：

bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed locus.bed

注意，locus.bed文件是必须的，示例如下：

chr1    10000   20000   gene1
chr1    20000   30000   gene2
chr1    30000   40000   gene3

bedtools getfasta

getfasta可以根据坐标区域来从基因组里面提取fasta序列

示例：

bedtools getfasta -fi genome.fa -bed HighQual_locus.bed -fo HighQual.fa

注意：bed格式用来记录坐标区域，参考基因组用-fi参数指定，输出的fasta序列文件用-fo参数指定

bedtools intersect

intersect可以分析特定位置上有哪些基因，接下来就还可以对基因进行各种下游分析。需要提供基因的坐标信息和需要处理的bed文件。

示例：

bedtools intersect -a features.bed -b locus.bed -wa -wb -F 1 | bedtools groupby -g 1-3 -c 7 -o collapse

注意：-f和-F可以用来指定覆盖比例(Minimum overlap required as a fraction of A/B)，详见下文或bedtools intersect -h。

常用bedtools intersect案例

用来求两个BED或者BAM文件中的overlap，overlap可以进行自定义是整个genome features的overlap还是局部。 加-wa参数可以报告出原始的在A文件中的feature，加-wb参数可以报告出原始的在B文件中的feature, 加-c参数可以报告出两个文件中的overlap的feature的数量，参数-s可以得到忽略strand的overlap，具体案例如下：

bedtools merge

merge用于合并位于同一个bed/gff/vcf文件中的重叠区域。

bedtools merge [OPTION] –i <bed/gff/vcf>

-s      必须相同(正负)链的区域才合并（软件默认不考虑正负链特征）
-n      报告合并的区域数量，没有被合并则1
-d      两个独立区域间距小于（等于）该值时将被合并为一个区域
-nms    报告被合并区域的名称
-scores 报告几个被合并特征区域的scores

其他小功能

参考资料

bcftools

BCFtools Manual

下载与安装

wget https://github.com/samtools/bcftools/releases/download/1.11/bcftools-1.11.tar.bz2
tar xjvf bcftools-1.11.tar.bz2
cd bcftools-1.11/
./configure --prefix=~/biosoft/bcftools
make
make install

bcftools annotate

annotate命令有两个用途:

bcftools concat

concat命令可以将多个VCF文件合并为一个VCF文件，要求输入的VCF文件必须是排序之后的，如果包含多个样本的信息，样本的顺序也必须一致。经典的应用场景包括合并不同染色体上的VCF文件，合并SNP和INDEL两种类型的VCF文件，用法如下:

bcftools concat a.vcf.gz b.vcf.gz -o concat.vcf

注意：输入文件必须是经过bgzip压缩的文件，而且还需要有tabix的索引。

bcftools merge

merge命令也是用于合并VCF文件，主要用于将单个样本的VCF文件合并成一个多个样本的VCF文件。用法如下:

bcftools merge a.vcf.gz b.vcf.gz -o merge.vcf

注意：

bcftools isec

isec用于在多个VCF文件之间取交集，差集，并集等操作，默认参数为取交集。经典的应用场景是对多种软件的SNP Calling结果进行Venn分析。用法如下, 其中-p参数指定输出结果的目录:

bcftools isec a.vcf.gz b.vcf.gz -p dir

bcftools stats

stats命令用于统计VCF文件的基本信息。如突变个数、突变类型的个数、转换颠换个数、测序深度、INDEL长度。还可以利用plot-vcfstats进行可视化处理。用法如下:

bcftools stats sample1.vcf > sample1.stats

输出文件可以用于plot-vcfstats命令进行可视化: plot-vcfstats sample1.stats -p output，这个脚本位于bcftools安装目录的misc目录下，依赖latex生成pdf文件。

bcftools index

index命令用于对VCF文件建立索引，要求输入的VCF文件必须是使用bgzip压缩之后的文件，支持.csi和.tbi两种索引，用法如下:

bcftools index sample1.vcf.gz       # csi格式
bcftools index -t sample1.vcf.gz    # tbi格式

bcftools view

bcftools query

query命令也是用于格式转换，和view命令不同，query通过表达式来指定输出格式，可定制化程度更高。用法如下:

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n' sample1.vcf.gz

-f参数通过一个表达式来指定输出格式，其中变量的写法如下:

bcftools sort

sort命令用于对VCF文件排序，按照染色体位置进行排序: bcftools sort sample1.vcf.gz -o sample1.sorted.vcf

bcftools reheader

reheader命令有两个用途:

常用例子：详见此文

BWA

BWA Manual
https://github.com/lh3/bwa

简介

BWA，即Burrows-Wheeler-Alignment Tool，是一种能够将差异度较小的序列比对到一个较大的参考基因组上的软件包。它有三个不同的算法：

对于上述三种算法，首先需要使用索引命令构建参考基因组的索引，用于后面的比对。所以，使用BWA整个比对过程主要分为两步，第一步建索引，第二步使用BWA MEM进行比对。

BWA的使用需要两中输入文件：

安装

git clone https://github.com/lh3/bwa.git
cd bwa
make

# quick start
./bwa index ref.fa
./bwa mem ref.fa read-se.fq.gz | gzip -3 > aln-se.sam.gz
./bwa mem ref.fa read1.fq read2.fq | gzip -3 > aln-pe.sam.gz

安装好之后，会出现一个名为bwa的可执行文件，输入下面命令可以查看帮助信息

./bwa

Program: bwa (alignment via Burrows-Wheeler transformation)
Version: 0.7.17-r1188
Contact: Heng Li <<a href="mailto:lh3@sanger.ac.uk" _href="mailto:lh3@sanger.ac.uk">lh3@sanger.ac.uk</a>>
Usage:   bwa <command> [options]
Command: index         index sequences in the FASTA format
         mem           BWA-MEM algorithm
         fastmap       identify super-maximal exact matches
         pemerge       merge overlapping paired ends (EXPERIMENTAL)
         aln           gapped/ungapped alignment
         samse         generate alignment (single ended)
         sampe         generate alignment (paired ended)
         bwasw         BWA-SW for long queries
         shm           manage indices in shared memory
         fa2pac        convert FASTA to PAC format
         pac2bwt       generate BWT from PAC
         pac2bwtgen    alternative algorithm for generating BWT
         bwtupdate     update .bwt to the new format
         bwt2sa        generate SA from BWT and Occ
Note: To use BWA, you need to first index the genome with `bwa index'.
      There are three alignment algorithms in BWA: `mem', `bwasw', and
      `aln/samse/sampe'. If you are not sure which to use, try `bwa mem'
      first. Please `man ./bwa.1' for the manual.

建立索引 index

bwa软件的作用是将 reads 比对到参考基因组上，在比对之前，首先需要对参考基因组建立索引，即对 fasta 文件构建 FM-index。

Usage:

bwa index [ –p prefix ] [ –a algoType ] <in.db.fasta>

# Example:
bwa index IRGSP-1.0_genome.fasta

Options:

BWA-MEM 算法

该算法先使用 MEM (maximal exact matches) 进行 seeding alignments，再使用 SW(affine-gap Smith-Waterman) 算法进行 seeds 的延伸。BWA–MEM 算法执行局部比对和剪接性。可能会出现 query 序列的多个不同的部位出现各自的最优匹配，导致 reads 有多个最佳匹配位点。有些软件如 Picard’s markDuplicates 跟 mem 的这种剪接性比对不兼容，在这种情况下，可以使用 –M 选项来将 shorter split hits 标记为次优。

Usage:

bwa mem [options] ref.fa reads_1.fq.gz [reads_2.fq.gz] > aln.sam

Options:

另：对于超长读长的reads,比如PacBio和Nanopore测序仪产生的reads, 用法如下:

bwa mem -x pacbio ref.fa reads.fq > aln.sam
bwa mem -x ont2d ref.fa reads.fq > aln.sam

BWA-backtrack 算法

对应的子命令为aln/samse/sample，先使用 aln 命令将单独的 reads 比对到参考序列，再使用 samse/sampe 生成 sam 文件。

Usage:

# single-end reads
bwa aln [options] ref.fa read.fq > aln_sa.sai
bwa samse [options] ref.fa aln_sa.sai read.fq > aln-se.sam

# paired-end reads
bwa aln [options] ref.fa read1.fq > aln_sa1.sai
bwa aln [options] ref.fa read2.fq > aln_sa2.sai
bwa sampe [options] ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln-pe.sam

Options:

BWA-SW 算法

Usage:

bwa bwasw -t <threads> ref.fa reads_1.fq.gz [reads_2.fq.gz] > aln.sam

有 2 个输入文件时，进行 paired-end 比对，此模式仅对 Illumina 的 short-insert 数据进行比对。在 Paired-end 模式下，BWA-SW 依然输出剪接性比对结果，但是这些结果会标记为 not properly paired; 同时如果有多个匹配位点，则不会写入 mate 的匹配位置。

aln mem bwasw怎么选

常用示例

BWA works with a variety types of DNA sequence data, though the optimal algorithm and setting may vary. The following list gives the recommended settings:

BWA-MEM is recommended for query sequences longer than ~70bp for a variety of error rates (or sequence divergence). Generally, BWA-MEM is more tolerant with errors given longer query sequences as the chance of missing all seeds is small. As is shown above, with non-default settings, BWA-MEM works with Oxford Nanopore reads with a sequencing error rate over 20%.

samtools

SAMtools Manual

samtools是一个用于操作sam和bam文件的工具合集。

samtools view

view命令的主要功能是：将输入文件转换成输出文件，通常是将比对后的sam文件转换为bam文件。；然后对bam文件进行各种操作，比如数据的排序(不属于本命令的功能)和提取(这些操作是对bam文件进行的，因而当输入为sam文件的时候，不能进行该操作)。

bam文件优点：bam文件为二进制文件，占用的磁盘空间比sam文本文件小；利用bam二进制文件的运算速度快。

view命令中，对sam文件头部的输入(-t或-T）和输出(-h)是单独的一些参数来控制的。

Usage:

samtools view [options] <in.bam>|<in.sam> [region1 [...]]

默认情况下不加 region，则是输出所有的 region

Options:

Examples:

（1）将sam文件和bam文件相互转换。注意：没有header的sam文件不能转换成bam文件。

samtools view -b -S lane.bam -o lane.sam
samtools view -h sample.bam > sample.sam

（2）提取比对到参考序列上的比对结果: -F

samtools view -bF 4 lane.bam > lane.F4.bam

（3）提取paired reads中两条reads都比对到参考序列上的比对结果，只需要把两个4+8的值12作为过滤参数即可

samtools view -bF 12 lane.bam > lane.F12.bam

（4）提取没有比对到参考序列上的比对结果: -f

samtools view -bf 4 lane.bam > lane.f4.bam

（5）提取paired reads中比对到参考基因组上的数据

samtools view -b -F 4 -f 8 lane.sam > only.read1.mapped.bam
samtools view -b -F 8 -f 4 lane.sam > only.read2.mapped.bam
samtools view -b -F 12 lane.sam > all.mapped.read1.read2.bam

（6）提取bam文件中比对到scaffold1上的比对结果，并保存到sam文件格式

samtools view lane.bam scaffold1 > scaffold1.sam

（7）提取scaffold1上能比对到30k到100k区域的比对结果

samtools view lane.bam scaffold1:30000-100000 > scaffold1_30k-100k.sam

（8）根据fasta文件，将 header 加入到 sam 或 bam 文件中

samtools view -T genome.fasta -h scaffold1.sam > scaffold1.h.sam

（9）sam 和 bam 格式转换

# BAM转换为SAM
samtools view -h 1.bam > out.sam

# SAM转换为BAM
samtools view -bS 1.sam > out.bam

# 如果无@SQ，要写-t，所以如果没有@SQ，需要：
samtools faidx ref.fa
samtools view -bt ref.fa.fai out.sam > out.bam

samtools sort

对bam文件进行排序。

Usage:

samtools sort [-l level] [-m maxMem] [-o out.bam] [-O format] [-n] [-T tmpprefix] [-@ threads] [in.sam|in.bam|in.cram]

Options:

samtools merge

当有多个样本的bam文件时，可以使用merge命令将这些bam文件进行合并为一个bam文件。merge命令将多个已经排序后的bam文件合并成为一个保持所有输入记录并保持现有排序顺序的bam文件。

Usage:

samtools merge [-nurlf] [-h inh.sam] [-R reg] [-b <list>] <out.bam> <in1.bam> [<in2.bam> ... <inN.bam>]

Options:

Examples:

#合并test_L1.bam和test_L2.bam文件
samtools merge -h test.sam \
    test_L1_L2.bam \
    test_L1.sorted.bam \
    test_L2.sroted.bam

#合并test_L1.bam和test_L2.bam文件中的指定区域chr7
samtools merge -h test.sam \
    -R chr7
    test_L1_L2.bam \
    test_L1.sorted.bam \
    test_L2.sroted.bam

samtools index

为了能够快速访问bam文件，可以为已经基于坐标排序后bam或者cram的文件创建索引，生成以.bai或者.crai为后缀的索引文件。必须使用排序后的文件，否则可能会报错。另外，不能对sam文件使用此命令。如果想对sam文件建立索引，那么可以使用tabix命令创建。

Usage:

samtools index [-bc] [-m INT] <aln.bam|cram> [out.index]

Options:

samtools faidx

对fasta文件建立索引，生成的索引文件以.fai后缀结尾。该命令也能依据索引文件快速提取fasta文件中的某一条（子）序列．如果没有指定区域，faidx命令就创建文件索引并生成后缀为.fai的索引文件。如果指定区域，便显示该区域序列。输入文件可以为bgzip压缩格式的文件。

Usage:

samtools faidx <ref.fa|ref.fa.gz>
samtools faidx <ref.fa|ref.fa.gz> <region1, ...> > selected.fa

生成的索引文件.fai是一个文本文件，共有5列：

samtools tview

tview能直观的显示出reads比对基因组的情况，和IGV有点类似。

Usage:

samtools tview <aln.bam> <ref.fa>

samtools flagstat <in.bam>

统计输入文件的相关数据并将这些数据输出至屏幕显示。每一项统计数据都由两部分组成，分别是QC pass和QC failed，表示通过QC的reads数据量和未通过QC的reads数量。以”PASS + FAILED”格式显示。

samtools depth

得到每个碱基位点或者区域的测序深度，并输出到标准输出。需要先index。

Usage:

samtools depth [options] [in1.bam|in1.sam|in1.cram[in2.bam|in2.sam|in2.cram]…]

Options:

samtools mpileup

该命令用于生成bcf文件，再使用bcftools进行SNP和Indel的分析。

Usage:

samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]

Options:

Examples:

samtools mpileup -f genome.fasta abc.bam > abc.txt
samtools mpileup -gSDf genome.fasta abc.bam > abc.bcf
samtools mpileup -guSDf genome.fasta abc.bam | bcftools view -cvNg - > abc.vcf
# complex one
samtools mpileup -l locus.bed -r 7 -q 1 -C 50 -t DP,DV -m 2 -F 0.002 -uvf human.fasta test.bam -o test.chr7.raw.vcf

mpileup不使用-u或-g参数时，则不生成二进制的bcf文件，而生成一个文本文件(输出到标准输出)。该文本文件统计了参考序列中每个碱基位点的比对情况。该文件每一行代表了参考序列中某一个碱基位点的比对结果，包含6列：参考序列名，位置，参考碱基，比对上的reads数，比对情况，比对上的碱基的质量。

其中第5列内容解释如下：

samtools rmdup

NGS上机测序前需要进行PCR一步，使一个模板扩增出一簇，从而在上机测序的时候表现出为1个点，即一个reads。若一个模板扩增出了多簇，结果得到了多个reads，这些reads的坐标(coordinates)是相近的。在进行了reads比对后需要将这些由PCR duplicates获得的reads去掉，并只保留最高比对质量的reads。使用rmdup命令即可完成。

Usage:

samtools rmdup [-sS] <input.sorted.bam> <output.bam>

Options:

samtools idxstats <in.bam>

检索和打印已建立索引的bam文件的统计数据，包括参考序列名称、序列长度、比对上的reads数量和未比对上的read数量。输出结果显示在屏幕上，以制表符分割。

samtools bedcov [-Q minQual] <region.bed> <in1.bam> [...]

计算由BED文件指定的基因组区域内的总碱基数量(read depth)。

samtools reheader <in.header.sam> <in.bam> > <out.bam>

替换bam文件的头

samtools cat [-h header.sam] [-o out.bam] <in1.bam> [...<inN.bam>]

连接多个bam文件，适用于非sorted的bam文件

将bam文件转换为fastq文件

有时候，我们需要提取出比对到一段参考序列的reads，进行小范围的分析，以利于debug等。这时需要将bam或sam文件转换为fastq格式。可以使用 bedtools bamtofastq -i <BAM> -fq <FQ1> -fq2 <FQ2> 完成。

GATK4

参考：
全基因组分析实践（参考：碱基矿工）
GATK4.0和全基因组数据分析实践（上）
GATK4.0和全基因组数据分析实践（下）

BWA_GATK_pipeline.sh

#!/bin/bash
# Basic BWA-GATK4 Pipeline for Resequencing (WGS) Analysis
# Author: liuyujie0136
# Date: 1/4/2022

## Install
conda install -y bwa samtools gatk4

## Go to Work Dir
cd ~/resequencing/
mv *.fq.gz original_fq/
mkdir renamed_fq/ sam_files/ bam_files/ gvcf_files/

## BWA Index
bwa index Ref.fa

## Rename FastQ Files
for i in `cat rename`
do
    ori=${i#*,}
    ren=${i%,*}
    
    ln -s ../original_fq/${ori}_R1.fq.gz renamed_fq/${ren}_R1.fq.gz
    ln -s ../original_fq/${ori}_R2.fq.gz renamed_fq/${ren}_R2.fq.gz
done

## BWA - mapping
for i in `cat rename`
do
    base=${i%,*}
    nohup bash 01.bwa_mem.sh $base Ref.fa renamed_fq/ sam_files/ "@RG\tID:${base}\tPL:illumina\tSM:${base}" &
done

## samtools - convert to bam and sort
for i in `cat rename`
do
    base=${i%,*}
    nohup bash 02.sam2bam.sh $base sam_files/ bam_files/ &
done

## GATK4 - mark duplicates and make GVCFs
gatk CreateSequenceDictionary -R Ref.fa

for i in `cat rename`
do
    base=${i%,*}
    nohup bash 03.gatk_make_gvcf.sh $base Ref.fa bam_files/ gvcf_files/ &
done

## GATK4 - Joint Calling of GVCFs and filter
nohup bash 04.gatk_call_variations.sh reseq_test Ref.fa gvcf_files/ ./ &

01.bwa_mem.sh

#!/bin/bash
# Usage: bash 01.bwa_mem.sh <FastQ Basename> <BWA Index> <Input FastQ Dir> <Output SAM Dir> <Read Group>
# Author: liuyujie0136
# Date: 1/4/2022

bwa mem -R "$5" -o $4/$1.sam $2 $3/${1}_R1.fq.gz $3/${1}_R2.fq.gz

Read Group 设置
e.g. @RG\tID:Sample1\tPL:illumina\tSM:Sample1

  1. ID，一般为测序的lane ID。
  2. PL，为测序的平台，需要准确填写。在GATK中，PL只被允许写为ILLUMINA,SLX,SOLEXA,SOLID,454,LS454,COMPLETE,PACBIO,IONTORRENT,CAPILLARY,HELICOS或UNKNOWN等信息之一。如果这一步没有设定正确的PL名称，则后续使用GATK过程则可能会出现报错。
  3. SM：为样本ID，因为我们在测序时，由于样本量较多，可能会分成许多不同的lane被分别测出来，这时候可以使用SM区分这些样本。
  4. LB：测序文库的名字。一般如果ID足够用于区分，则可以不设定LB。
    对于序列比对而言，设定这四个参数就足够了。除此之外，在RG中还需要用制表符\t将各个内容区分开。

02.sam2bam.sh

#!/bin/bash
# Usage: bash 02.sam2bam.sh <SAM Basename> <Input SAM Dir> <Output BAM Dir>
# Author: liuyujie0136
# Date: 1/4/2022

samtools view -b -h -o $3/$1.bam $2/$1.sam

samtools sort -o $3/$1.sorted.bam $3/$1.bam

rm $3/$1.bam

samtools index $3/$1.sorted.bam

03.gatk_make_gvcf.sh

#!/bin/bash
# Usage: bash 03.gatk_make_gvcf.sh <BAM Basename> <Ref File> <Input BAM Dir> <Output GVCF Dir>
# Author: liuyujie0136
# Date: 1/4/2022

gatk MarkDuplicates -I $3/$1.sorted.bam -O $3/$1.sorted.markdup.bam -M $3/$1.sorted.markdup.metrics

samtools index $3/$1.sorted.markdup.bam

gatk HaplotypeCaller -R $2 -I $3/$1.sorted.markdup.bam -O $4/$1.g.vcf.gz -ERC GVCF

04.gatk_call_variations.sh

#!/bin/bash
# Usage: bash 04.gatk_call_variations.sh <VCF Basename> <Ref File> <Input GVCF Dir> <Output VCF Dir>
# Author: liuyujie0136
# Date: 1/4/2022

gvcf=""
for i in `ls $3/*.g.vcf.gz`
do
    gvcf=$gvcf"-V $i "
done

gatk CombineGVCFs -R $2 $gvcf -O $3/$1.combined.g.vcf.gz

gatk GenotypeGVCFs -R $2 -V $3/$1.combined.g.vcf.gz -O $4/$1.vcf.gz

gatk VariantFiltration -V $4/$1.vcf.gz \
    --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
    --filter-name "LowConf" \
    -O $4/$1.filtered.vcf.gz

bowtie与bowtie2

Bowtie Manual
Bowtie2 Manual

vcftools

VCFtools Manual

vcftools群体遗传参数计算示例

1. Tajima’s D，这个是选择相关的一个参数，大于0代表群体观测杂合度高于预期杂合度，稀有等位基因频率降低（群体收缩或者平衡选择），小于0说明群体观测杂合位点少于预期值，稀有等位基因频率增加（群体扩张或者低频选择）。 也就是说，只有0是正常的，其他都是选择发生。

vcftools \
  --vcf test.vcf \
  --TajimaD 500000 \
  --out taj

2. π，核苷酸多样性，越大说明核苷酸多样性越高，越低说明两个座位DNA序列差异越小。

vcftools \
  --vcf test.vcf \
  --window-pi 500000 \
  --out pi

记录一个关于计算核酸多样性（pi）的经历
https://blog.csdn.net/weixin_43362619/article/details/107239012
之前在做叶绿体基因组核酸多样的时候，先是用全基因组做，选择窗口和步长使用DnaSP5来求pi值。后来发现文章里放的都是编码序列和非编码序列，或者每个基因的pi值。叶绿体的编码序列一般在80个多个，如果再加上非编码区，那么有一百多个序列需要计算，手动算的话挺费时费劲的，一分钟算三个的话得近一个小时…作为懒人的我肯定是不愿意每天花一个小时做这样的重复工作，于是…
我找了文章里的计算方法，发现使用的竟然全部是DnaSP，后来发现vcftools好像可以计算pi，而且是linux版本的，通过脚本很容易实现批量计算，于是抓紧试了试。但是vcftools只支持vcf格式的文件，没关系，把比对的序列转成vcf格式就行了，最终实现批量计算多个基因的pi值，小小开心下。但是有个问题，同样的序列用vcftools计算得到的pi和用DnaSP计算得到的pi不一样。开始是觉得我转换得到vcf的时候不对。于是我找到了另一款计算pi的工具：perl中有个Bio::PopGen::Statistics模块，可以计算pi，而且输入的是比对后的序列。但是！结果仍然和DnaSP计算的不一样！令人惊讶的是，popgene和vcftools的计算结果是相同的，我百思不得其解。难道是DnaSP算的有问题。于是找到了计算pi的公式，手动算了下，结果和DnaSP的值相同。难道…另外两个算错了？于是研究了Bio::PopGen::Statistics模块中计算pi值的代码，终于找到了原因。
Bio::PopGen::Statistics在计算pi的时候和vcftools类似，先把比对结果转换成vcf格式，然后计算。但是不同的是，如果两个序列比如AAA,AAT，直接计算的话差异碱基数为1。当转成vcf的时候是用 3 A T 0/0 1/1，来表示(3号位点其中一个是A，另一个是T)。这时候的差异碱基数就变成了4（把0/0，1/1分别看成两个位点，如果还原成序列的话就是A，A，T，T，及从原来的2条序列变成了4条…），实际计算的时候是把样品数当成4而不是2。也就是说，差异碱基数变成了原先的四倍，样品（序列）数变成了原先的两倍，所以最终算的值和DnaSP的结果不同。
考虑到vcftools是用于vcf文件的计算，如果有个位点的信息是0/1，那么确实要表示成两个碱基，所以样品数量要变成原先的两倍，计算的也没毛病。DnaSP是直接计算序列的pi，没有考虑太多，和直接套公式计算的结果相同，结果没得说。Bio::PopGen::Statistics呢，我给的是比对的序列，然后你按照vcf的来算，似乎并不友好（没详细看它的文档），按理说我给它的序列和给DnaSP的序列是相同的，结果应该也是相同的，但是实际情况并不是。
综上所述，如果想计算的是比对后序列的pi，DnaSP首选，在这里也是唯一的选择，因为即使你把序列给Bio::PopGen::Statistics算，它算的结果依旧不是你想要的…如果你的结果是群体的变异位点信息，那么就选择vcftools吧~（当然，不怕麻烦的话也可以把变异位点还原成序列，然后用DnaSP来做…）。

3. Fst, 分化系数，从0到1说明亲缘关系越来越远。接近于0说明两个个体亲缘关系近，接近1说明亲缘关系远。10number.txt和20number.txt是想要分组的个体文件，每行一个个体，跟vcf中的个体名称一致。

vcftools \
  --vcf test.vcf \
  --weir-fst-pop pop1.txt \
  --weir-fst-pop pop2.txt \
  --out pop1vspop2 \
  --fst-window-size 500000

4. Hardy-Weinberg平衡检测，这个主要是检测基因型频率是否等于基因频率乘积。比如A:0.3,a:0.7那么Aa的频率是否为0.42

vcftools \
  --vcf test.vcf \
  --hardy \
  --out hw

5. XP-CLR

https://github.com/hardingnj/xpclr

这个是计算群体分化的一个指标，主要是基于选择清除（selective sweep）信息进行群体分化的计算。个人感觉比Fst更新的地方是用了多个SNP对应的等位基因频率估计染色体一段区域的分化情况，而Fst是用单个SNP对应的等位基因频率。相对来说，XP-CLR估计群体分化会更准确一些。

xpclr \
  --format vcf \
  --input test.vcf \
  --samplesA pop1.txt \
  --samplesB pop2.txt \
  --chr 1 \
  --minsnps 20 \
  --maxsnps 600 \
  --size 10000 \
  --step 10000 \
  --out out.chr1.xpclr