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
注意要点：
- -sn =sample name 样本名（选择目标样本）
- -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列 （不一致按突变类型筛选）
选择“I列 ExonicFunc” 筛选非同义突变
选择“Data”
选择“Filter”
选择“autoFilter”
去除“0”
去除“unknown”
去除“synonymous SNV”
得到非同义突变等突变有意义的位点
新建“sheet3”
把“sheet2”中全选，
拷贝到“sheet3”中
把“sheet2”中“I列”全选
确认
选择“sheet2”中“Z列scsnv” 这一列是可剪切的，在内含子里，前面突变类型怕漏或在最后一个外显子编码区或同义突变
选择“Data”
选择“Filter”
选择“autoFilter”
去除“0”
得到剪切位点有意义的位点
把“sheet2”中全选，
拷贝到“sheet3”中原有数据中下面 得到候选位点了
选择“sheet3”中“K列 genomic super” 根据基因表型关联得到基因这个表型是如何的
选择“Data”
选择“Filter”
选择“autoFilter”
选择“0”
得到没有同源序列的基因
“G”列就是用于Phenolzer的genelist
第二种方法
把WJ-2302-candidate.hg19_multianno.txt转成WJ-2338-candidate.hg19_multianno.xlsx
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
网友评论