WES流程圖.PNG

測序數(shù)據(jù)的輸入文件為raw fastq，質(zhì)控軟件為fastqc 和trim_galore，通過參數(shù)設(shè)置，得到結(jié)果文件clean fastq。

比對。輸入文件為clean/HQ fastq，比對軟件使用BWA（不止這一種），比對算法MEM，結(jié)果文件為raw sam，然后轉(zhuǎn)成bam文件。

找變異。IGV可視化bam，GATK過濾bam，GATK找變異得到vcf

注釋。輸入文件為HQ vcf，比對軟件Annovar，使用參考數(shù)據(jù)庫，得到結(jié)果文件即與功能相關(guān)的變異信息。

后續(xù)進行個性化分析。

實戰(zhàn)階段

軟件安裝

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

更改環(huán)境變量

echo 'export PATH=/home/jianmingzeng/biosoft/myBin/bin:$PATH' >>~/.bashrc 
這個/home/jianmingzeng/biosoft/myBin/bin是個例子，得要找到自己安裝了的那個軟件的位置。怎么找到？以安裝GTK軟件為例，先 ./gatk --help 能出來，然后 pwd 出來路徑，例如 /home/vip18/GATK4/gatk-4.0.11.0  然后將這個添加到環(huán)境變量里。最后這么使之生效：
source ~/.bashrc

這事也可以這么干：
進入vim編輯，
vim ~/.bashrc
寫進去，
export PATH="/home/jianmingzeng/biosoft/myBin/bin:$PATH"
使之生效，
source ~/.bashrc

配置鏡像

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

創(chuàng)建名為wes的軟件安裝環(huán)境

conda create -n wes python=2 #這里同時安裝了python 2版本的軟件。也可以后續(xù)安裝它。

查看當前conda環(huán)境

conda info --envs

激活/進入conda的wes環(huán)境（避免每次用-n wes）

source activate wes

安裝sra-tools軟件

conda search sra-tools   # fastq-dump --help
conda install -y sra-tools   # done正確安裝，且能調(diào)出軟件help

繼續(xù)安裝軟件

conda install -y fastqc trim-galore star bwa hisat2 bowtie2 subread tophat htseq bedtools deeptools salmon cutadapt multiqc # 注意核查每個軟件是否能夠被調(diào)用

數(shù)據(jù)下載

找到SRA數(shù)據(jù)

https://www.ncbi.nlm.nih.gov/sra?term=SRP077971 [可修改SRP號]

找到SRA數(shù)據(jù)集，找到_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

source activate wes # step1 激活環(huán)境

cat SRR_Acc_List.txt | while read id; do (prefetch ${id});done  # step2獲得SRR_Acc_List.txt 或 觀察數(shù)據(jù)編號的規(guī)律下載

ls ./1.SRA/*sra | while read id; do ( nohup fastq-dump --gzip --split-3 -O ./ ${id} ); done # step3 SRA數(shù)據(jù)轉(zhuǎn)成fastq

fastq 格式

fastq格式.PNG

現(xiàn)在用預(yù)先下載好的文件做練習(xí)，例如
7E5239.L1_1.fastq.gz
7E5241.L1_2.fastq.gz

ll ./*/*fastq.gz  # 可以搜索找到gz結(jié)尾的文件。

fastqc軟件對練習(xí)數(shù)據(jù)進行質(zhì)控

fastqc -t 3 ./7E5239.L1_1.fastq.gz  ./7E5241.L1_2.fastq.gz

質(zhì)控結(jié)果會反映在每一個.html文件里，樣品較多看起來不方便，可使用multiqc整合在一個.html文件里。multiqc不止能整合fastqc軟件運行結(jié)果，也可以整合其他軟件運行結(jié)果。更多信息需官網(wǎng)查詢。

運行fastqc時候可能會報錯，例如出現(xiàn)Permission denied 字樣，表明沒有相應(yīng)權(quán)限，這說明該軟件結(jié)果的輸出不在自己的目錄中。查詢幫助即 fastqc --help 查找出入?yún)?shù)用法，將輸出路徑更改到自己的目錄中來，即定義輸出路徑。例如：
fastqc /home/qmcui/test_boy/1.fastq_qc/*.fastq.gz -o ./project/wes/

所有樣品的fastqc運行成功后再運行multiqc。可先multiqc --help，查看幫助文檔。如下運行：
multiqc ./

查看結(jié)果，得知reads大小，作為參考，以備運行trim_galore過濾數(shù)據(jù)，實際該步驟通常不需要自己做，過程較為復(fù)雜，通常是測序公司做好即可，以下trim_galore的運行僅為了解。

trim_galore --phred33 -q 25 -e 0.1--length 36 --stringency 3--paired -o ./ 
./7E5241.L1_1.fastq.gz ./7E5241.L1_2.fastq.gz

批量質(zhì)控。多個樣本時候需要批量質(zhì)控，寫腳本文件后綴為.sh，例如取名為trim_galore.sh

cat > 自己的目錄 ./trim_galore.sh   # 例如，$ cat /home/qmcui/test_boy/1.fastq_qc/trim_galore.sh 

# 以下為寫入文本的內(nèi)容
cat config|while read id
do 
    arr=(${id})
    fq1=${arr[0]}
    fq2=${arr[1]}
    echo "################################################" # 設(shè)定這個是為調(diào)試程序。另外，這一長串比較顯眼，容易在屏幕上找到。這并不重要。
    trim_galore -q 25 --phred33 --length 36 \
    --stringency 3 --paired -o ./  $fq1 $fq2
    echo "OK"
done

# crtl +c結(jié)束

# 在運行之前，先cat 一下剛才寫的文件，再查看一下所寫內(nèi)容正確與否，確認無誤后再運行。

# 提交后臺運行
nohup bash trim_galore.sh & v# 會返回一個后臺任務(wù)id號

jobs -l   #可以看到剛才返回的后臺任務(wù)id已經(jīng)在運行

Mapping

比對目的：
? 將打斷測序的reads比回參考基因組
? 得到比對結(jié)果sam文本，用于后續(xù)分析

比對策略：
? index先建立索引
? mem算法(maximal exact matches)

軟件：
? bwa(Burrows-Wheeler Aligner )

參考基因組：
? hg38.fa

Sort

sort目的：
? 用于后續(xù)軟件分析
? 提取bam序列

比對策略：

一步到位生成sort.bam （bam文件排序）
? bwa生成sam轉(zhuǎn)bam,然后再排序。
? bwa比對時直接利用管道符‘|’直接生成排序bam

軟件：
? samtools(sort參數(shù))

一步到位生成sort.bam

bwa mem -t 5 -R "@RG\tlD:$sample\tSM:$sample\tLB:WGS\tPL:illumina" $INDEX $fq1 $fa2 |samtools sort -@ 5 -o $sample.sort.bam -
 

bam統(tǒng)計

ls *.bam|while read id;do(samtools flagstat $id > $(basename$id".bam").stat);done

Mark PCR重復(fù)(Mark Duplicates)

Mark Duplicates目的：
? 標記/刪除PCR重復(fù)的reads
? 為后續(xù)call變異位點增加可信度,去掉假陽性

使用軟件GATK4

建立路徑

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

找變異——variant call之gatk多樣本

參數(shù)：
? sample賦值
? -I 輸入文件
? -O 輸出文件
? 1> 重定向標準輸出
? 2> 重定向錯誤輸出到1的文件內(nèi)(目的：記錄程序運行日志)

目的：
? 得到vcf前矯正的bam步驟

軟件：
? GATK 子命令FixMateInformation
? samtools 子命令index 建索引

step 1

### 建立路徑 ###
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

step 2

目的：
? 得到vcf前矯正的bam步驟

參數(shù)：
? sample賦值樣本ID
? ref賦值參考基因組
? snp賦值dbsnp數(shù)據(jù)庫(同indel)
? -I 輸入文本
? -O 輸出文本(承前接后)

軟件：
? 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

step 3

參數(shù)：
? sample賦值樣本ID
? ref賦值參考基因組
? snp賦值dbsnp數(shù)據(jù)庫
? -I 輸入文本
? -O 輸出文本(承前接后)
目的：
? 得到vcf前矯正的bam步驟

軟件——GATK 子命令
? ApplyBQSR
? HaplotypeCaller

### 矯正bam，先賦值 ####
sample="7E5241"
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

step 4

### calling variation得到raw vcf ###
mkdir gatk_single_vcf && cd gatk_single_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

raw vcf filter

步驟較多，暫時跳過

注釋 （使用annovar）

# 建立路徑（沒有的話，要先建一個）
# mkdir annovar

~qmcui/software/annovar/annovar/convert2annovar.pl \
-format vcf4old \      # -filter  # -allsample  #另有許多參數(shù)，詳見官網(wǎng)
${sample}_filtered.vcf > ${sample}.annovar

# 運行完成之后得到.anno.variant_function的文件，可以對這個文件進行進一步的操作，例如：

cat  7E5241.anno.variant_function

less -S  7E5241.anno.variant_function|cut -f 2 | sort|uniq -c #  統(tǒng)計第二列

less -S  7E5241.anno.variant_function|grep -v 'unknown'|grep -v -w 'synonymous SNV'|cut -f 2 | sort|uniq -c

less -S  7E5241.anno.variant_function|grep -v 'unknown'|grep -v -w 'synonymous SNV'|less -S >##.txt

less -S 7E5240_last_annovar.exonic_variant_function|grep -v 'unknown'|grep -v -w 'synonymous SNV'|awk '$4=="chr1" && $5 >0 && $6 <208559422 {print $0}'|less -S
`