下载参考基因组
wget https://ftp.ncbi.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.28_GRCh38.p13/GCA_000001405.28_GRCh38.p13_genomic.fna.gz
1. gzip -d GCA_000001405.28_GRCh38.p13_genomic.fna.gz
下载gff注释文件
wget https://ftp.ncbi.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.28_GRCh38.p13/GCA_000001405.28_GRCh38.p13_genomic.gff.gz

构建索引

samtools faidx GCA_000001405.28_GRCh38.p13_genomic.fna
bwa index GCA_000001405.28_GRCh38.p13_genomic.fna
java -jar picard.jar CreateSequenceDictionary R=GCA_000001405.28_GRCh38.p13_genomic.fna

比对

bwa mem
-t 2 \ # -t 线程数
-R '@RG\tID:J2\tSM:J2\tPL:illumina' \ #添加read group信息
/home/ug0983/peng/GCA_000001405.28_GRCh38.p13_genomic.fna \ # 参考基因组index路径
/home/ug0983/peng/Sample_JZ23031529-HY0625-FFPE-2-J2/JZ23031529-HY0625-FFPE-2-J2_combined_R1.fastq.gz \ # fq1数据路径
/home/ug0983/peng/Sample_JZ23031529-HY0625-FFPE-2-J2/JZ23031529-HY0625-FFPE-2-J2_combined_R2.fastq.gz \ # fq2数据路径
2>J2.bwa.log | \ # 比对日志
samtools
sort -@ 2 \ # 线程数
-m 1G \ # 每个线程内存大小
-o S1.sort.bam - # -o 指定输出文件名称

bwa mem -t 2 -R '@RG\tID:J2\tSM:J2\tPL:illumina' /home/ug0983/peng/GCA_000001405.28_GRCh38.p13_genomic.fna /home/ug0983/peng/Sample_JZ23031529-HY0625-FFPE-2-J2/JZ23031529-HY0625-FFPE-2-J2_combined_R1.fastq.gz /home/ug0983/peng/Sample_JZ23031529-HY0625-FFPE-2-J2/JZ23031529-HY0625-FFPE-2-J2_combined_R2.fastq.gz
2> J2.bwa.log | samtools sort -@ 2 -m 1G -o S1.sort.bam -
samtools view S1.sort.bam|less -S

去重复

java
-Xmx4g -XX:ParallelGCThreads=2 -jar /home/ug0983/pengt/picard.jar
MarkDuplicates I=S1.sort.bam \ # 输入文件名称
O=S1.sort.markdup.bam \ # 输出文件名称
CREATE_INDEX=true \ # 创建bam index文件
REMOVE_DUPLICATES=true \ # dup reads是否删除
M=S1.marked_dup_metrics.txt

java -Xmx4g -XX:ParallelGCThreads=2 -jar /home/ug0983/peng/picard.jar MarkDuplicates I=S1.sort.bam O=S1.sort.markdup.bam CREATE_INDEX=true REMOVE_DUPLICATES=true M=S1.marked_dup_metrics.txt

比对率统计

samtools flagstat S1.sort.bam > S1.sort.bam.flagstat
cat S1.sort.bam.flagstat

覆盖度统计

samtools coverage S1.sort.bam > S1.sort.bam.coverage
cat S1.sort.bam.coverage

多组样品比对

准备好样品名称及测序数据路径文件

$ cat sample.list（文档需要再liunix中创建，windows下的文档不兼容））
S1 /home/ug0983/peng/J24/J2_combined_R1.fastq.gz /home/ug0983/peng/J24/J2_combined_R2.fastq.gz
S2 /home/ug0983/peng/J24/J4_combined_R1.fastq.gz /home/ug0983/peng/J24/J4_combined_R2.fastq.gz

创建批量比对脚本（文档需要再liunix中创建，windows下的文档不兼容）

!/bin/bash

ref=/home/ug0983/peng/GCA_000001405.28_GRCh38.p13_genomic.fna
cat sample.list |while read sp fq1 fq2
do
echo "
bwa mem -t 2 -R '@RG\tID: $sp\tSM:$ sp\tPL:illumina' $ref$ fq1 $fq2 2> J2.bwa.log | samtools sort -@ 2 -m 1G -o$ sp.sort.bam -

java -Xmx4g -XX:ParallelGCThreads=2 -jar /home/ug0983/peng/picard.jar MarkDuplicates I= $sp.sort.bam O=$ sp.sort.markdup.bam CREATE_INDEX=true REMOVE_DUPLICATES=true M=$sp.marked_dup_metrics.txt

samtools flagstat $sp.sort.bam >$ sp.sort.bam.flagstat

samtools coverage $sp.sort.bam >$ sp.sort.bam.coverage

" > mapping.$sp.sh
done

人类参考基因组下载比对