1. 原始數(shù)據(jù)（fastq） → 質(zhì)量控制（fastqc）

操作及路徑：
在fastqc路徑下雙擊打開(kāi)可以修改腳本餐抢，設(shè)置參數(shù)剩愧，右鍵open in terminal后運(yùn)行腳本：./fastqc
/home/wei/Documents/gatk/FastQC/fastqc
點(diǎn)擊File→ open→ Documents→ align
注意要點(diǎn)：
1）數(shù)據(jù)量
數(shù)據(jù)量=reads數(shù)目*reads的長(zhǎng)度肮之；
1G的數(shù)據(jù)量=十億個(gè)堿基
2）每個(gè)reads質(zhì)量分布(除了測(cè)序深度注意數(shù)據(jù)均一性即覆蓋度，遺傳咨詢分析胚系突變，通常20x>98%合格)

2. 比對(duì)基因組（speedseq）bam文件

路徑：點(diǎn)擊File→open→Documents→align
雙擊打開(kāi)align.sh，修改腳本參數(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后運(yùn)行腳本：
sh -v align.sh
注意要點(diǎn)：
1）每次運(yùn)行前修改align.sh腳本里與樣本有關(guān)的信息如tSM:k373
2）align為比對(duì)命令（樣本測(cè)序數(shù)據(jù)與hg19數(shù)據(jù)庫(kù)進(jìn)行比對(duì)）
3）-t 線程數(shù)（線程越多，數(shù)據(jù)分析速度越快）
4）-o 輸出結(jié)果文件名稱
5）-T 數(shù)據(jù)結(jié)果存放路徑（一般為當(dāng)前路徑殴泰，也可寫絕對(duì)路徑）
6）-R 將比對(duì)的reads進(jìn)行分組
7）該分析流程只針對(duì)雙端測(cè)序數(shù)據(jù)于宙，不支持單端測(cè)序數(shù)據(jù)（k373-1.fastq.gz k373-2.fastq.gz原始數(shù)據(jù)）

3. 計(jì)算覆蓋度與深度

路徑：點(diǎn)擊Documents→day1→coverage.sh
腳本內(nèi)容如下：
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后運(yùn)行腳本：
sh -v coverage.sh
用excel打開(kāi)查看數(shù)據(jù)質(zhì)量結(jié)果：Documents→day1→test1.sample_summary
目標(biāo)區(qū)域數(shù)據(jù)量/總的數(shù)據(jù)量=中靶率
遺傳樣本數(shù)據(jù)：20x以上的區(qū)域占目標(biāo)區(qū)域百分之多少， >98%為合格
注意要點(diǎn)：
1）-Xmx3g中的3g為跑該程序所允許的最大內(nèi)存悍汛，根據(jù)自己電腦性能修改
2）-T 測(cè)序深度算法命令調(diào)用
3）-R 參考基因組數(shù)據(jù)庫(kù)
4）-o 輸出結(jié)果文件名稱
5）-I 輸入數(shù)據(jù)捞魁，待計(jì)算測(cè)序深度的樣本文件
6）-L 要分析的基因組位置、目標(biāo)區(qū)域（捕獲區(qū)域：指定比對(duì)的數(shù)據(jù)庫(kù)离咐，如exome.bed全外范圍內(nèi)谱俭，exome-chr1.bed為1號(hào)染色體全外范圍內(nèi)）
7）-ct 測(cè)序深度x占比（10x，20x宵蛀，50x昆著，根據(jù)遺傳、腫瘤等樣本需求可自己設(shè)置）

4. 檢測(cè)變異（GATK）

待分析樣本與參考基因組相比發(fā)生突變的术陶，得到vcf（這個(gè)位點(diǎn)變異成什么堿基凑懂，基因型是雜合型還是純合型，reads怎么樣梧宫，有多少條reads支持該位點(diǎn)的突變接谨，位點(diǎn)突變峰度，突變深度塘匣，位點(diǎn)深度多高脓豪，突變類型是否準(zhǔn)確）
分析思路：一次可分析多個(gè)樣本變異信息（包括單個(gè)或多個(gè)先證者、家系對(duì)照組等）→ 低質(zhì)量突變位點(diǎn)去除→ 待分析樣本突變數(shù)據(jù)從總數(shù)據(jù)中提取→ 反向提取對(duì)照樣本數(shù)據(jù)（用來(lái)參考的未突變數(shù)據(jù)）→ 待分析樣本的獨(dú)有突變位點(diǎn)（對(duì)照樣本擁有的相同位點(diǎn)會(huì)被過(guò)濾）
路徑：點(diǎn)擊Documents→day1→call.sh
修改腳本參數(shù)設(shè)置
運(yùn)行sh -v call.sh

call.sh腳本內(nèi)容如下：

1.一次分析多個(gè)樣本變異信息

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
注意要點(diǎn)：
1） -Xmx6g中的6g表示用最大6G內(nèi)存去跑程序
2） -T UnifiedGenotyper 分析樣本中哪些位置與參考基因組不一致
3） -I 輸入bam文件（list里可把要比對(duì)的bam都放進(jìn)去）
4） -L 捕獲區(qū)域（指定比對(duì)的數(shù)據(jù)庫(kù)范圍忌卤，如exome.bed全外顯子范圍內(nèi)扫夜，exome-chr1.bed為1號(hào)染色體全外顯子范圍內(nèi)）
5） -glm BOTH 檢測(cè)SNP（點(diǎn)突變）和INDEL（微小插入缺失50bp內(nèi)插入或重復(fù)缺失）
6） -stand_emit_conf 10 高于10分的位置才檢測(cè)變異
-stand_call_conf 30 高于30分的位置認(rèn)為可信
分?jǐn)?shù)影響因素：1）測(cè)序深度 2）突變的reads占總reads的比率、 每條reads比對(duì)得分
7）打開(kāi)test1.vcf→vcf 表頭含有正文里短語(yǔ)縮寫的解釋驰徊，運(yùn)行的命令历谍，參考基因組大小，染色體有多長(zhǎng)等信息辣垒。正文：染色體編號(hào)望侈，位置信息，參考基因組堿基勋桶，突變堿基脱衙，變異得分侥猬，變異的可靠性評(píng)估
8）GT:AD:DP:GQ:PL 1/1:0,101:101:99:3328,247,0
GT對(duì)應(yīng)突變基因型：1/1 純合突變，0/1雜合突變 捐韩，0/0野生型退唠，./.該位置沒(méi)有reads覆蓋
AD對(duì)應(yīng)野生型和突變型覆蓋數(shù)目：0,101即檢測(cè)到0個(gè)野生型reads，檢測(cè)到101次突變型reads（真正雜合子數(shù)目差不多荤胁，0.5）
DP對(duì)應(yīng)測(cè)序深度瞧预，101次（多少條reads覆蓋該位點(diǎn)）
GQ對(duì)應(yīng)基因型得分，如99仅政，分?jǐn)?shù)越高表示該突變檢測(cè)數(shù)據(jù)越可靠
PL該虛擬機(jī)內(nèi)部評(píng)估模型參數(shù)垢油，與測(cè)序數(shù)據(jù)質(zhì)量相關(guān)。
9）變異得分影響因素：測(cè)序深度圆丹，突變r(jià)eads占總reads的比率滩愁，每條reads比對(duì)得分測(cè)序深度太低會(huì)出現(xiàn)假陽(yáng)性
10）Documents→day1→test1.vcf.idx文件為test1.vcf的縮影文件

2.去除樣本突變信息中低質(zhì)量位點(diǎn)

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
注意要點(diǎn)：
1）-T SelectVariants 選擇變異位點(diǎn)，GATK篩選位點(diǎn)的子命令辫封。
2）-ef 去除低質(zhì)量的變異位點(diǎn)（參數(shù)設(shè)置默認(rèn)低于30分的位點(diǎn)會(huì)被過(guò)濾）

3.待分析樣本突變數(shù)據(jù)從總數(shù)據(jù)中提取

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
注意要點(diǎn)：

1. -sn =sample name 樣本名（選擇目標(biāo)樣本）
2. -env 去除該樣本中的野生型的位點(diǎn)（只保留有突變的位點(diǎn)硝枉，只提取有變異的，野生型位點(diǎn)會(huì)被過(guò)濾）

4.反向提取對(duì)照樣本（對(duì)照樣本不能為相同疾病樣本倦微，未判斷明確遺傳模式的情況下也不能有血緣關(guān)系）

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
注意要點(diǎn)：
1）-xl_sn 反向選擇除了指定樣本之外的樣本
2）若對(duì)照樣本為相同疾病樣本妻味，有可能將該類疾病共同突變位點(diǎn)過(guò)濾掉，造成假陰性
3）未判斷明確遺傳模式的情況下欣福，有血緣關(guān)系樣本弧可，可能會(huì)將該家族共同隱性突變位點(diǎn)過(guò)濾，造成假陰性
4）隨機(jī)選擇正常健康人數(shù)據(jù)作為對(duì)照樣本

5.待分析樣本的獨(dú)有突變位點(diǎn)

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
注意要點(diǎn)：
1）--discordance 選擇與其他參考樣本vcf（WJ-2338-other.vcf）不一致的位點(diǎn)
2）對(duì)照樣本中任何一個(gè)樣本含有與待分析樣本相同的位點(diǎn)劣欢，該位點(diǎn)都會(huì)被去掉棕诵，縮小篩選位點(diǎn)的數(shù)目
3）對(duì)照數(shù)不建議設(shè)太多，直接刪會(huì)導(dǎo)致假陰性結(jié)果凿将，刪3-4個(gè)最常見(jiàn)snv去掉校套，同源性的假陽(yáng)性位點(diǎn)去掉

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
注意要點(diǎn)：
1）-buildver vcf所用的基因組版本
2）-remove 去除注釋過(guò)程中的中間文件
3）-protocol 使用哪些數(shù)據(jù)庫(kù)作為對(duì)照來(lái)注釋
4）-nastring . 用點(diǎn)號(hào)替代缺省的值，該部分?jǐn)?shù)據(jù)為參考數(shù)據(jù)庫(kù)中均未收錄牧抵，屬于研究未報(bào)道位點(diǎn)
5）-vcfinput 注釋文件格式是vcf
g,r,r,f,f,f,f,f,f,f,f,f,f,f 各個(gè)數(shù)據(jù)庫(kù)用什么樣的方式向這些位點(diǎn)在數(shù)據(jù)庫(kù)上進(jìn)行關(guān)聯(lián)
6）-operation 用什么方式來(lái)注釋笛匙，注釋結(jié)果輸出格式（g,r,f, , ,）
g 表示輸出結(jié)果按照gene匹配，輸出
r 表示按基因堿基片段注釋（突變和什么疾病相關(guān)犀变，疾病和基因相關(guān)妹孙，基因在基因組上有一段區(qū)間范圍，
如chr1染色體第10000-20000的一段區(qū)域获枝，只要在這區(qū)域里任何一個(gè)位點(diǎn)突變蠢正，直接歸到這個(gè)疾病）這個(gè)點(diǎn)在這個(gè)區(qū)域里就把它注釋成某個(gè)疾病省店，點(diǎn)的區(qū)域進(jìn)行關(guān)聯(lián)
f 點(diǎn)對(duì)點(diǎn)注釋（按堿基匹配嚣崭，千人基因組頻率）

6.篩選位點(diǎn)（Excel）

過(guò)濾思路（遺傳病的致病位點(diǎn)）：通過(guò)人群頻率1000g2014oct_all和exac_eas小于0.005（千分之5）→看clinvar數(shù)據(jù)庫(kù)有沒(méi)有已經(jīng)報(bào)道的致病的或可能致病的位點(diǎn)笨触，與患者表型符合的，如有雹舀，寫報(bào)告評(píng)估位點(diǎn)致病性芦劣，如沒(méi)有，根據(jù)突變類型進(jìn)行篩選非同義突變说榆，去除“0”虚吟、去除“unknown”、去除“synonymous SNV”（這樣留下來(lái)的非同義突變签财、錯(cuò)位突變串慰、移碼突變、插入與缺失導(dǎo)致的非移碼突變荠卷，這些突變類型才能影響基因的功能）→選scsnv一列模庐，判斷是否影響它的剪切的烛愧，不為0油宜，（會(huì)影響它的剪切）→形成候選reads

那么多候選位點(diǎn)哪一類基因與患者的表型符合
通過(guò)Phenolzer（在線或者本地）將基因和表型進(jìn)行關(guān)聯(lián)
在線：用bing搜索Phenolzer，輸入郵箱怜姿，輸入患者表型（Hpo）或者疾采髟（Als）圈越大顏色越深代表和這個(gè)表型越相關(guān)
本地：用篩選后的基因列表（disease輸入表型，candidate-gene-list輸入候選基因列表沧卢，pheno腳本運(yùn)行進(jìn)行關(guān)聯(lián)）→還能通過(guò)XYAutoFilterV4進(jìn)行篩選→看候選基因列表漸凍人有哪些基因與漸凍人相關(guān)蚁堤，用到panelapp這個(gè)網(wǎng)站（非常系統(tǒng)整理了目前所有單基因遺傳病的基因列表、表型等）→把基因列表下的基因名稱拷貝到一個(gè)txt文本里→用XYAutoFilterV4篩選在這個(gè)疾病下罕見(jiàn)非同義突變

WJ-2338-candidate.hg19_multianno文件用excel打開(kāi)但狭，另存為xlsx格式
A列Chr 染色體編號(hào)
B列Start 變異起始的位置
C列End 終止的位置
D列Ref 參考基因組的堿基
E列Alt 突變了什么堿基
F列Func.refGene 變異所處參考基因的功能區(qū)（外顯子編碼區(qū)披诗，UTR區(qū)，內(nèi)含子區(qū)立磁，基因鍵區(qū)）
G列Gene.refGene 變異所處參考基因名稱（如果是基因間呈队，則是兩側(cè)的基因）
H列GeneDetail.refGene 非翻譯區(qū) （到距離有多遠(yuǎn)）*174G>C,
終止密碼子往后174位G>C→ 3'翻譯區(qū)終止密碼子在后面
如果是5'，就要往前C-
I列ExonicFunc.refGene 突變類型（ 同義突變唱歧，非移碼的插入宪摧，錯(cuò)位突變）
nonsynonymous SNV（非同義突變）在這里只能是錯(cuò)位突變
J列AAChange.refGene 氨基酸水平的轉(zhuǎn)變，轉(zhuǎn)錄轉(zhuǎn)移颅崩，幾號(hào)外顯子几于，編碼區(qū)從哪里到哪里，（同一個(gè)基因可能具有多個(gè)轉(zhuǎn)錄本沿后，氨基酸改變的位置在不同的轉(zhuǎn)錄本中有可能不一樣）
k列g(shù)enomicSuperDups 位于基因的同源區(qū) 不為0（很有可能假陽(yáng)性位點(diǎn)）
L列phastConsElements46way 保守性 不為0為保守沿彭，多個(gè)物種中都是保守的，序列上保守尖滚，這個(gè)地方如突變更可能會(huì)致病
M列esp6500siv2_all 6500個(gè)人群頻率
N列ExAC_Freq 65000個(gè)人群頻率
O-U列ExAC_AFR膝蜈、ExAC_AMR锅移、ExAC_EAS、ExAC_FIN饱搏、ExAC_NFE非剃、ExAC_OTH、ExAC_SAS 65000個(gè)人分別各個(gè)人種 用最多的EAS東亞（6000人左右）
V列1000g2014oct_eas 千人基因組東亞 （500人左右）
W列 1000g2014oct_all 所有的千人基因組2500個(gè)人
ExAC和1000g區(qū)別
1000g是全基因組推沸，有很多內(nèi)含子突變备绽，ExAC外顯子水平，1000g全基因組水平鬓催，內(nèi)含子都有
X列avsnp142 比對(duì)位點(diǎn)有個(gè)rs號(hào)
Y列clinvar_20140929 Clinvar數(shù)據(jù)庫(kù)被收錄就會(huì)在這里注釋不同級(jí)別 （良性還是致病的肺素、可能致病的、這種解讀有差異的宇驾、解讀有沖突等）
clinvar_dba 和什么疾病相關(guān)
Z列scSNV 可變剪切數(shù)據(jù)庫(kù) （有些位點(diǎn)外顯子和內(nèi)含子交接的地方突變倍靡，會(huì)影響剪切，不為0可能會(huì)影響剪切）
AA课舍、AB列revel mcapv1.0 最新版錯(cuò)位突變塌西，預(yù)測(cè)錯(cuò)位突變危害性的2個(gè)數(shù)據(jù)庫(kù)（不為0代表有害）
AC列cosmic68wgs 是腫瘤的
AD列至最后列SIFT_score、SIFT_pred 各種算法對(duì)非同義突變保守性預(yù)測(cè)值(預(yù)測(cè)錯(cuò)位突變) T：耐受性筝尾；D：有害捡需，Polyphen2_HDIV_pred（D：可能有害，P：可能有害筹淫；B：良性）
最后是vcf文件基因型

A列染色體號(hào)全選右鍵點(diǎn)擊Column Width...
width選擇1
選擇“Q”列到“W”列
選擇“DATA”
選擇“FILTER”
選擇“Standred Filter”
選擇“1000g2014oct_all” and ”<“ and”0.005”
and
選擇“exac_eas” and ”<“ and”0.005”
點(diǎn)擊“+”增加一個(gè)”sheet2”
把”test1-clean.hg19.multianno”內(nèi)容全選后
拷到“sheet2”中

選擇“Y列 clinvar2014” 看報(bào)道的和病人表型是符合的
選擇“Data”
選擇“Filter”
選擇“autoFilter”
選擇“pathogenic” 數(shù)據(jù)庫(kù)報(bào)道的致病或可能致病的站辉， 挑出來(lái)還把這些位點(diǎn)它的表型和你bam的表型是否一致的
看這個(gè)clinvar描述的這個(gè)基因關(guān)聯(lián)的疾病是否與這個(gè)樣本的表型一致
如果一致，就找到了結(jié)果
如果不一致损姜，重新全選Y列 （不一致按突變類型篩選）

把“sheet2”中“I列”全選
確認(rèn)
選擇“sheet2”中“Z列scsnv” 這一列是可剪切的，在內(nèi)含子里摧阅，前面突變類型怕漏或在最后一個(gè)外顯子編碼區(qū)或同義突變
選擇“Data”
選擇“Filter”
選擇“autoFilter”
去除“0”
得到剪切位點(diǎn)有意義的位點(diǎn)
把“sheet2”中全選汰蓉，
拷貝到“sheet3”中原有數(shù)據(jù)中下面 得到候選位點(diǎn)了

python XYAutoFilterV4.py

(Genelist = 'als.txt' #可選項(xiàng)：.txt文件的genelist,把als.txt刪了）

練習(xí) (****篩選漸凍人的基因）
通過(guò)panelsapp這個(gè)網(wǎng)站去尋找和漸凍人相關(guān)基因的列表，這些基因下載到本地逸尖，復(fù)制粘貼到genelist文本里古沥，提交XYAutoFilterV4，只分析漸凍人相關(guān)基因上罕見(jiàn)位點(diǎn)

打開(kāi)Documents→vcf→anno （2338漸凍人）

新建文件夾 test 把WJ-2338-candidate.hg19_multianno.xlsx和XYAutoFilterV4.py粘貼進(jìn)去

打開(kāi)XYAutoFilterV4.py 娇跟，Genelist = 'als.txt' #可選項(xiàng)：.txt文件的genelist,把a(bǔ)ls.txt刪了岩齿，save.

運(yùn)行python XYAutoFilterV4.py

打開(kāi)wj2338-resultals3.xlsx→把H列Gene.refGene下的基因復(fù)制
去網(wǎng)頁(yè)搜索bing→phenolyzer（http://phenolyzer.wglab.org/）→Diseases/Phenotypes（輸入疾病全稱）（如als漸凍人 Amyotrophic lateral sclerosis） →Gene Selection選yes→Enter your genes here （輸入 wj2338-resultals.xlsx的H列）→點(diǎn)擊submit→點(diǎn)擊網(wǎng)址自動(dòng)跳出

第二種方法：

打開(kāi)Documents→phenolyzer

填寫 disease →疾病全稱（如als漸凍人 Amyotrophic lateral sclerosis）

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

sh -v pheno.sh

結(jié)果在out文件夾里

打開(kāi)網(wǎng)址 https://panelapp.genomicsengland.co.uk/

點(diǎn)擊 panels

輸入 疾病全稱（如als漸凍人 Amyotrophic lateral sclerosis）

下載

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

打開(kāi) XYAutoFilterV4.py

Genelist = 'als.txt' #可選項(xiàng)：.txt文件的genelist,把a(bǔ)ls.txt加上

python XYAutoFilterV4.py