- 下载参考基因组
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- 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\tPL:illumina'
fq1
sp.sort.bam -
java -Xmx4g -XX:ParallelGCThreads=2 -jar /home/ug0983/peng/picard.jar MarkDuplicates I=sp.sort.markdup.bam CREATE_INDEX=true REMOVE_DUPLICATES=true M=$sp.marked_dup_metrics.txt
samtools flagstat sp.sort.bam.flagstat
samtools coverage sp.sort.bam.coverage
" > mapping.$sp.sh
done
