1. 原始数据（fastq） → 质量控制（fastqc）

操作及路径：
在fastqc路径下双击打开可以修改脚本，设置参数，右键open in terminal后运行脚本：./fastqc
/home/wei/Documents/gatk/FastQC/fastqc
点击File→ open→ Documents→ align
注意要点：
1）数据量
数据量=reads数目*reads的长度；
1G的数据量=十亿个碱基
2）每个reads质量分布(除了测序深度注意数据均一性即覆盖度，遗传咨询分析胚系突变，通常20x>98%合格)

2. 比对基因组（speedseq）bam文件

路径：点击File→open→Documents→align
双击打开align.sh，修改脚本参数：
/home/wei/bin/speedseq align -t 1 -o k373 -R "@RG\tID:id\tSM:k373\tLB:lib" -T 'pwd' /home/wei/Documents/gatk/library/ucsc.hg19.fasta k373-1.fastq.gz k373-2.fastq.gz
右键open in terminal后运行脚本：
sh -v align.sh
注意要点：
1）每次运行前修改align.sh脚本里与样本有关的信息如tSM:k373
2）align为比对命令（样本测序数据与hg19数据库进行比对）
3）-t 线程数（线程越多，数据分析速度越快）
4）-o 输出结果文件名称
5）-T 数据结果存放路径（一般为当前路径，也可写绝对路径）
6）-R 将比对的reads进行分组
7）该分析流程只针对双端测序数据，不支持单端测序数据（k373-1.fastq.gz k373-2.fastq.gz原始数据）

3. 计算覆盖度与深度

路径：点击Documents→day1→coverage.sh
脚本内容如下：
java -Xmx3g -jar /home/wei/bin/GenomeAnalysisTK.jar -T DepthOfCoverage -R /home/wei/Documents/gatk/library/ucsc.hg19.fasta -o test1 -I WJ-2338.bam -L /home/wei/Documents/gatk/library/Agilent-v6.bed --omitDepthOutputAtEachBase --omitIntervalStatistics -ct 1 -ct 10 -ct 20 -ct 50
右键open in terminal后运行脚本：
sh -v coverage.sh
用excel打开查看数据质量结果：Documents→day1→test1.sample_summary
目标区域数据量/总的数据量=中靶率
遗传样本数据：20x以上的区域占目标区域百分之多少， >98%为合格
注意要点：
1）-Xmx3g中的3g为跑该程序所允许的最大内存，根据自己电脑性能修改
2）-T 测序深度算法命令调用
3）-R 参考基因组数据库
4）-o 输出结果文件名称
5）-I 输入数据，待计算测序深度的样本文件
6）-L 要分析的基因组位置、目标区域（捕获区域：指定比对的数据库，如exome.bed全外范围内，exome-chr1.bed为1号染色体全外范围内）
7）-ct 测序深度x占比（10x，20x，50x，根据遗传、肿瘤等样本需求可自己设置）

4. 检测变异（GATK）

待分析样本与参考基因组相比发生突变的，得到vcf（这个位点变异成什么碱基，基因型是杂合型还是纯合型，reads怎么样，有多少条reads支持该位点的突变，位点突变峰度，突变深度，位点深度多高，突变类型是否准确）
分析思路：一次可分析多个样本变异信息（包括单个或多个先证者、家系对照组等）→ 低质量突变位点去除→ 待分析样本突变数据从总数据中提取→ 反向提取对照样本数据（用来参考的未突变数据）→ 待分析样本的独有突变位点（对照样本拥有的相同位点会被过滤）
路径：点击Documents→day1→call.sh
修改脚本参数设置
运行sh -v call.sh

call.sh脚本内容如下：

1.一次分析多个样本变异信息

java -jar -Xmx6g /home/wei/bin/GenomeAnalysisTK.jar -T UnifiedGenotyper -R /home/wei/Documents/gatk/library/ucsc.hg19.fasta -I bam.list -L /home/wei/Documents/gatk/library/exome-chr1.bed -glm BOTH -stand_emit_conf 10 -stand_call_conf 30 -o test1.vcf
注意要点：
1） -Xmx6g中的6g表示用最大6G内存去跑程序
2） -T UnifiedGenotyper 分析样本中哪些位置与参考基因组不一致
3） -I 输入bam文件（list里可把要比对的bam都放进去）
4） -L 捕获区域（指定比对的数据库范围，如exome.bed全外显子范围内，exome-chr1.bed为1号染色体全外显子范围内）
5） -glm BOTH 检测SNP（点突变）和INDEL（微小插入缺失50bp内插入或重复缺失）
6） -stand_emit_conf 10 高于10分的位置才检测变异
-stand_call_conf 30 高于30分的位置认为可信
分数影响因素：1）测序深度 2）突变的reads占总reads的比率、 每条reads比对得分
7）打开test1.vcf→vcf 表头含有正文里短语缩写的解释，运行的命令，参考基因组大小，染色体有多长等信息。正文：染色体编号，位置信息，参考基因组碱基，突变碱基，变异得分，变异的可靠性评估
8）GT:AD:DP:GQ:PL 1/1:0,101:101:99:3328,247,0
GT对应突变基因型：1/1 纯合突变，0/1杂合突变 ，0/0野生型，./.该位置没有reads覆盖
AD对应野生型和突变型覆盖数目：0,101即检测到0个野生型reads，检测到101次突变型reads（真正杂合子数目差不多，0.5）
DP对应测序深度，101次（多少条reads覆盖该位点）
GQ对应基因型得分，如99，分数越高表示该突变检测数据越可靠
PL该虚拟机内部评估模型参数，与测序数据质量相关。
9）变异得分影响因素：测序深度，突变reads占总reads的比率，每条reads比对得分测序深度太低会出现假阳性
10）Documents→day1→test1.vcf.idx文件为test1.vcf的缩影文件

2.去除样本突变信息中低质量位点

java -Xmx6g -jar /home/wei/bin/GenomeAnalysisTK.jar -R /home/wei/Documents/gatk/library/ucsc.hg19.fasta -T SelectVariants -V test1.vcf -ef -o test1-clean.vcf
注意要点：
1）-T SelectVariants 选择变异位点，GATK筛选位点的子命令。
2）-ef 去除低质量的变异位点（参数设置默认低于30分的位点会被过滤）

3.待分析样本突变数据从总数据中提取

java -Xmx5g -jar /home/wei/bin/GenomeAnalysisTK.jar -R /home/wei/Documents/gatk/library/ucsc.hg19.fasta -T SelectVariants -V test1-clean.vcf -sn WJ-2338 -env -o WJ-2338.vcf
注意要点：

1. -sn =sample name 样本名（选择目标样本）
2. -env 去除该样本中的野生型的位点（只保留有突变的位点，只提取有变异的，野生型位点会被过滤）

4.反向提取对照样本（对照样本不能为相同疾病样本，未判断明确遗传模式的情况下也不能有血缘关系）

java -Xmx5g -jar /home/wei/bin/GenomeAnalysisTK.jar -R /home/wei/Documents/gatk/library/ucsc.hg19.fasta -T SelectVariants -V test1-clean.vcf -xl_sn WJ-2338 -env -o WJ-2338-other.vcf
注意要点：
1）-xl_sn 反向选择除了指定样本之外的样本
2）若对照样本为相同疾病样本，有可能将该类疾病共同突变位点过滤掉，造成假阴性
3）未判断明确遗传模式的情况下，有血缘关系样本，可能会将该家族共同隐性突变位点过滤，造成假阴性
4）随机选择正常健康人数据作为对照样本

5.待分析样本的独有突变位点

java -Xmx5g -jar /home/wei/bin/GenomeAnalysisTK.jar -R /home/wei/Documents/gatk/library/ucsc.hg19.fasta -T SelectVariants -V WJ-2338.vcf --discordance WJ-2338-other.vcf -env -o WJ-2338-candidate.vcf
注意要点：
1）--discordance 选择与其他参考样本vcf（WJ-2338-other.vcf）不一致的位点
2）对照样本中任何一个样本含有与待分析样本相同的位点，该位点都会被去掉，缩小筛选位点的数目
3）对照数不建议设太多，直接删会导致假阴性结果，删3-4个最常见snv去掉，同源性的假阳性位点去掉

5.变异注释（Annovar）

mkdir anno
/home/wei/annovar2/table_annovar.pl WJ-2338-candidate.vcf /home/wei/annovar2/humandb/hg19/ -buildver hg19 -out anno/WJ-2338-candidate -remove -protocol refGene,genomicSuperDups,phastConsElements46way,esp6500siv2_all,exac03,1000g2014oct_eas,1000g2014oct_all,avsnp142,clinvar_20180603,scsnv,revel,mcap,cosmic68wgs,ljb26_all -operation g,r,r,f,f,f,f,f,f,f,f,f,f,f -nastring . -vcfinput
注意要点：
1）-buildver vcf所用的基因组版本
2）-remove 去除注释过程中的中间文件
3）-protocol 使用哪些数据库作为对照来注释
4）-nastring . 用点号替代缺省的值，该部分数据为参考数据库中均未收录，属于研究未报道位点
5）-vcfinput 注释文件格式是vcf
g,r,r,f,f,f,f,f,f,f,f,f,f,f 各个数据库用什么样的方式向这些位点在数据库上进行关联
6）-operation 用什么方式来注释，注释结果输出格式（g,r,f, , ,）
g 表示输出结果按照gene匹配，输出
r 表示按基因碱基片段注释（突变和什么疾病相关，疾病和基因相关，基因在基因组上有一段区间范围，
如chr1染色体第10000-20000的一段区域，只要在这区域里任何一个位点突变，直接归到这个疾病）这个点在这个区域里就把它注释成某个疾病，点的区域进行关联
f 点对点注释（按碱基匹配，千人基因组频率）

6.筛选位点（Excel）

过滤思路（遗传病的致病位点）：通过人群频率1000g2014oct_all和exac_eas小于0.005（千分之5）→看clinvar数据库有没有已经报道的致病的或可能致病的位点，与患者表型符合的，如有，写报告评估位点致病性，如没有，根据突变类型进行筛选非同义突变，去除“0”、去除“unknown”、去除“synonymous SNV”（这样留下来的非同义突变、错位突变、移码突变、插入与缺失导致的非移码突变，这些突变类型才能影响基因的功能）→选scsnv一列，判断是否影响它的剪切的，不为0，（会影响它的剪切）→形成候选reads

那么多候选位点哪一类基因与患者的表型符合
通过Phenolzer（在线或者本地）将基因和表型进行关联
在线：用bing搜索Phenolzer，输入邮箱，输入患者表型（Hpo）或者疾病（Als）圈越大颜色越深代表和这个表型越相关
本地：用筛选后的基因列表（disease输入表型，candidate-gene-list输入候选基因列表，pheno脚本运行进行关联）→还能通过XYAutoFilterV4进行筛选→看候选基因列表渐冻人有哪些基因与渐冻人相关，用到panelapp这个网站（非常系统整理了目前所有单基因遗传病的基因列表、表型等）→把基因列表下的基因名称拷贝到一个txt文本里→用XYAutoFilterV4筛选在这个疾病下罕见非同义突变

WJ-2338-candidate.hg19_multianno文件用excel打开，另存为xlsx格式
A列Chr 染色体编号
B列Start 变异起始的位置
C列End 终止的位置
D列Ref 参考基因组的碱基
E列Alt 突变了什么碱基
F列Func.refGene 变异所处参考基因的功能区（外显子编码区，UTR区，内含子区，基因键区）
G列Gene.refGene 变异所处参考基因名称（如果是基因间，则是两侧的基因）
H列GeneDetail.refGene 非翻译区 （到距离有多远）*174G>C,
终止密码子往后174位G>C→ 3'翻译区终止密码子在后面
如果是5'，就要往前C-
I列ExonicFunc.refGene 突变类型（ 同义突变，非移码的插入，错位突变）
nonsynonymous SNV（非同义突变）在这里只能是错位突变
J列AAChange.refGene 氨基酸水平的转变，转录转移，几号外显子，编码区从哪里到哪里，（同一个基因可能具有多个转录本，氨基酸改变的位置在不同的转录本中有可能不一样）
k列genomicSuperDups 位于基因的同源区 不为0（很有可能假阳性位点）
L列phastConsElements46way 保守性 不为0为保守，多个物种中都是保守的，序列上保守，这个地方如突变更可能会致病
M列esp6500siv2_all 6500个人群频率
N列ExAC_Freq 65000个人群频率
O-U列ExAC_AFR、ExAC_AMR、ExAC_EAS、ExAC_FIN、ExAC_NFE、ExAC_OTH、ExAC_SAS 65000个人分别各个人种 用最多的EAS东亚（6000人左右）
V列1000g2014oct_eas 千人基因组东亚 （500人左右）
W列 1000g2014oct_all 所有的千人基因组2500个人
ExAC和1000g区别
1000g是全基因组，有很多内含子突变，ExAC外显子水平，1000g全基因组水平，内含子都有
X列avsnp142 比对位点有个rs号
Y列clinvar_20140929 Clinvar数据库被收录就会在这里注释不同级别 （良性还是致病的、可能致病的、这种解读有差异的、解读有冲突等）
clinvar_dba 和什么疾病相关
Z列scSNV 可变剪切数据库 （有些位点外显子和内含子交接的地方突变，会影响剪切，不为0可能会影响剪切）
AA、AB列revel mcapv1.0 最新版错位突变，预测错位突变危害性的2个数据库（不为0代表有害）
AC列cosmic68wgs 是肿瘤的
AD列至最后列SIFT_score、SIFT_pred 各种算法对非同义突变保守性预测值(预测错位突变) T：耐受性；D：有害，Polyphen2_HDIV_pred（D：可能有害，P：可能有害；B：良性）
最后是vcf文件基因型

A列染色体号全选右键点击Column Width...
width选择1
选择“Q”列到“W”列
选择“DATA”
选择“FILTER”
选择“Standred Filter”
选择“1000g2014oct_all” and ”<“ and”0.005”
and
选择“exac_eas” and ”<“ and”0.005”
点击“+”增加一个”sheet2”
把”test1-clean.hg19.multianno”内容全选后
拷到“sheet2”中

选择“Y列 clinvar2014” 看报道的和病人表型是符合的
选择“Data”
选择“Filter”
选择“autoFilter”
选择“pathogenic” 数据库报道的致病或可能致病的， 挑出来还把这些位点它的表型和你bam的表型是否一致的
看这个clinvar描述的这个基因关联的疾病是否与这个样本的表型一致
如果一致，就找到了结果
如果不一致，重新全选Y列 （不一致按突变类型筛选）

把“sheet2”中“I列”全选
确认
选择“sheet2”中“Z列scsnv” 这一列是可剪切的，在内含子里，前面突变类型怕漏或在最后一个外显子编码区或同义突变
选择“Data”
选择“Filter”
选择“autoFilter”
去除“0”
得到剪切位点有意义的位点
把“sheet2”中全选，
拷贝到“sheet3”中原有数据中下面 得到候选位点了

python XYAutoFilterV4.py

(Genelist = 'als.txt' #可选项：.txt文件的genelist,把als.txt删了）

练习 (****筛选渐冻人的基因）
通过panelsapp这个网站去寻找和渐冻人相关基因的列表，这些基因下载到本地，复制粘贴到genelist文本里，提交XYAutoFilterV4，只分析渐冻人相关基因上罕见位点

打开Documents→vcf→anno （2338渐冻人）

新建文件夹 test 把WJ-2338-candidate.hg19_multianno.xlsx和XYAutoFilterV4.py粘贴进去

打开XYAutoFilterV4.py ，Genelist = 'als.txt' #可选项：.txt文件的genelist,把als.txt删了，save.

运行python XYAutoFilterV4.py

打开wj2338-resultals3.xlsx→把H列Gene.refGene下的基因复制
去网页搜索bing→phenolyzer（http://phenolyzer.wglab.org/）→Diseases/Phenotypes（输入疾病全称）（如als渐冻人 Amyotrophic lateral sclerosis） →Gene Selection选yes→Enter your genes here （输入 wj2338-resultals.xlsx的H列）→点击submit→点击网址自动跳出

第二种方法：

打开Documents→phenolyzer

填写 disease →疾病全称（如als渐冻人 Amyotrophic lateral sclerosis）

candidate-gene-list →wj2338-resultals.xlsx的H列

sh -v pheno.sh

结果在out文件夹里

打开网址 https://panelapp.genomicsengland.co.uk/

点击 panels

输入 疾病全称（如als渐冻人 Amyotrophic lateral sclerosis）

下载

Amyotrophic lateral sclerosis_motor neuron disease.tsv修改名字为 als.txt

打开 XYAutoFilterV4.py

Genelist = 'als.txt' #可选项：.txt文件的genelist,把als.txt加上

python XYAutoFilterV4.py