wes实战流程
1.软件安装
- miniconda安装
wget -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda2-4.5.11-Linux-x86_64.sh
bash Miniconda2-4.5.11-Linux-x86_64.sh
- 配置镜像
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes
- 创建名为wes的软件安装环境
conda create -n wes python=2
- 查看当前conda环境
conda into --envs
- 激活/进入conda的wes环境，避免每次用-n wes
source activate wes
- 安装 sra-tools软件
conda search sra-tools #fastq-dump --help
conda install -y sra-tools #done正确安装，且能调成软件help
.....
- 退出环境
source deactivate #退出当前环境
2.数据下载
- 找到SRA数据
https://www.ncbi.nlm.nih.gov/sra?term=SRP077971 [可修改SRP号]
- 找到SRA数据集，找到_List.txt
https://www.ncbi.nlm.nih.gov/Traces/study/?WebEnv=NCID_1_5569658_130.14.22.76_5555_1540039 270_2975166503_0MetA0_S_HStore&query_key=1
-
获得fastq
1.激活环境
source activate wes
2.获得SRR_Acc_List.txt或观察数据编号的规律下载
cat SRR_Acc_List.txt | while read id; do (prefetch ${id});done
3. SRA数据转成fastq
ls ./1.SRA/*sra | while read id; do ( nohup fastq-dump --gzip --split-3 -O ./ ${id} ); done
Fastq格式
1.id 2.sequende/reads 3.+ 4.quanlity score
1.质控
mkdir 1.fastq_qc #./当前路径
cd 1.fastq_qc
###1.data statistics
fastqc -t 3 ./7E5241.L1_1.fastq.gz \
./7E5241.L1_2.fastq.gz# 可继续加*fq,有些软件不能用通配符
multiqc ./ # 整合结果
###2.filter data
trim_galore --phred33 \
-q25 -e0.1\
--length 36 --stringency 3 \
--paired -o ./ \
./7E5241.L1_1.fastq.gz \
./7E5241.L1_2.fastq.gz
批量指控
# qc.sh脚本的内容
cat $1 |while read id
do arr=(${id})
fq1=${arr[0]}
fq2=${arr[1]}
echo $fq1 $fq2
$bin_trim_galore -q 25 --phred33 \
--length 36 --stringency 3 --paired \
-o $dir $fq1 $fq2
done
# shell若报错解决报错
# source deactivate
2.比对Mapping
Mapping
1.比对目的：1）将打断测序的reads比回参考基因组
2）得到比对结果sam文本，用于后续分析
2.比对策略：1）index先建立索引
2）mem算法（maximal exact matches）
3.软件：bwa（Burrows-Wheeler Aligner）
4.参考基因组：hg38.fa
###建立project路径
mkdir 0.config 2.mapping
###1.make index[可使用别人建好的]
cd 0.config
ln -s /public/biosoft/GATK/resources/bundle/hg38/Homo_sa piens_assembly38.fasta `./hg38.fa`
bwa index ./hg38.fa
###2.single sample mapping(nohup ... &)
cd 2.mapping
bwa mem -t 5 \
-R "@RG\tID:7E5241\tSM:7E5241\tLB:WGS\tPL:Illumina" \ ../0.config/hg38 \ ../1.fastq_qc/result_dir/7E5241.L1_1_val_1.fq.gz \ ../1.fastq_qc/result_dir/7E5241.L1_2_val_2.fq.gz \ 1>7E5241.sam 2>7E5241.log
bwa结果
sam格式
1."\t"分割每列
2.@是头文件
3.比对行=必需11列+一列可变字段
4.与bam内容一模一样
bam特点
1.bam文件为其二进制文件
2.站内存小
3.一般存bam删sam
- sam转bam
samtool view -bS -h \
7E5241.sam -o 7E5241.bam
# 或:
samtool view -bS -h 7E5241.sam > 7E5241.bam
- bam查看方法
samtool view -h smaple.bam |less -S 2 samtool view smaple.bam |less -S
3.bam文件排序-sort
sort
1.sort目的：1）用于后续软件分析
2）提取bam序列
2.比对策略： 1）bwa生产sam转bam，然后再排序
2）bwa比对时直接利用管道符“|”直接生成排序bam
3.软件：samtools（sort参数）
madir 3.sort_bam #建立路径
cd 3.sort_bam
sort bam
samtools sort -@ 5 ../2.mapping/7E5241.bam \
-o 7E5241.sort.bam
###一步到位生成sort.bam
bwa mem -t 5 -R
"@RG\tID:$sample\tSM:$sample\tLB:WGS\tPL:Illumina" $INDEX $fq1 $fq2 | samtools sort -@ 5 -o $sample.sort.bam -
###bam统计
ls *.bam | while read id ;do (samtools flagstat $id > $(basename $id ".bam").stat);done
4.Mark PCR重复
1.Mark Duplicates目的：1）标记/删除PCR重复的reads
2）为后续call变异位点增加可信度，去掉假阳性
2.软件：GATK(MarkDuplicates)
mkdir 4.gatk_markdup #建立路径
cd 4.gatk_markdup
mkdup bam donot delete
sample="7E5241"
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" MarkDuplicates \
-I ../3.sort_bam/$sample.sort.bam \ -O ${sample}_marked.bam \
-M ${sample}.metrics \
1>${sample}_log.mark 2>&1
5.找变异-variant call之gatk多样本
1.目的：得到vcf前矫正的bam步骤
2.参数：• sample赋值
• -I 输入文件
• -O 输出文件
• 1> 重定向标准输出
• 2> 重定向错误输出到1的文件内（目的：记录程序运行日志）
3.软件：1）GATK子命令FixMateInformation
2）samtools子命令index建索引
mkdir 6.gatk_bam #建立路径
cd 6.gatk_bam
sample="7E5241"
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" FixMateInformation \
-I ${sample}_marked.bam \
-O ${sample}_marked_fixed.bam \ -SO coordinate \
1>${sample}_log.fix 2>&1
samtools index ${sample}_marked_fixed.bam
1.目的：• 得到vcf前矫正的bam步骤
2.参数：• sample赋值样本ID
• ref赋值参考基因组
• snp赋值dbsnp数据库(同indel)
• -I 输入文本
• -O 输出文本(承前接后)
3.软件：• GATK 子命令BaseRecalibrator
cd 6.gatk_bam # 进入路径
sample="7E5241" # 找变异，先赋值ref="/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly 38.fasta" snp="/public/biosoft/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz" indel="/public/biosoft/GATK/resources/bundle/hg38/Mills_and_1000G_gold _standard.indels.hg38.vcf.gz"
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" BaseRecalibrator \ -R $ref \
-I ${sample}_marked_fixed.bam \
--known-sites $snp \
--known-sites $indel \
-O ${sample}_recal.table
1.目的：• 得到vcf前矫正的bam步骤
2.参数：• sample赋值样本ID
• ref赋值参考基因组
• snp赋值dbsnp数据库
• -I 输入文本
• -O 输出文本(承前接后)
3.软件：• GATK 子命令ApplyBQSR
• HaplotypeCaller
sample="7E5241" #矫正bam，先赋值
ref="/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly 38.fasta"
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" ApplyBQSR
-R $ref \
-I ${sample}_marked_fixed.bam \ -bqsr ${sample}_recal.table \
-O ${sample}_bqsr.bam \
1>${sample}_log.ApplyBQSR 2>&1
###calling variation 得到vcf
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" HaplotypeCaller \
-R $ref \
-I ${sample}_bqsr.bam \
--dbsnp $snp \
-O ${sample}_raw.vcf \
1>${sample}_log.HC 2>&1
image
image
6注释
annovar
1.目的：注释vcf
2.参数：1）convert2annovar perl脚本生成注释需要文件
2）annotata_variation注释
3.软件：annovar的perl脚本
###建立路径
mkdir annovar
cd annovar
###注释
sample="7E5241"
~qmcui/software/annovar/annovar/convert2annovar.pl \
-format vcf4old \
${sample}_filtered.vcf > ${sample}.annovar
~qmcui/software/annovar/annovar/annotate_variation.pl \
-buildver hg38 \
--outfile ${sample}.anno \
${sample}.annovar \ /public/biosoft/ANNOVAR/annovar/humandb
image
网友评论