vcftools --gzvcf input.vcf --chr n --recode – recode-INFO-all --stdout | gzip -c output.vcf.gz
我們提供的服務(wù)有：成都做網(wǎng)站、網(wǎng)站建設(shè)、微信公眾號開發(fā)、網(wǎng)站優(yōu)化、網(wǎng)站認(rèn)證、門頭溝ssl等。為數(shù)千家企事業(yè)單位解決了網(wǎng)站和推廣的問題。提供周到的售前咨詢和貼心的售后服務(wù)，是有科學(xué)管理、有技術(shù)的門頭溝網(wǎng)站制作公司
說明：
–gzvcf：處理壓縮格式的vcf文件（可替換為–vcf）
–chr n：選擇染色體n，例：–chr 1
–recode：重新編碼為vcf文件,有過濾操作都要加上--recode
–recode-INFO-all：將輸出的文件保存所有INFO信息
–stdout：標(biāo)準(zhǔn)輸出，后接管道命令
–gzip -c：壓縮
--max-missing
--max-missing的取值是0-1，為1時表示某個位點(diǎn)上所有的樣本必須都有基因型，一個樣本的基因型都不能缺。所以這個選項(xiàng)可以理解為：能分型的樣本占總樣本的比例至少為多少。
基本的思想就是利用數(shù)據(jù)流重定向，把原來輸出到屏幕上的數(shù)據(jù)定向""到文件里
ATAC-seq信息分析流程主要分為以下幾個部分：數(shù)據(jù)質(zhì)控、序列比對、峰檢測、motif分析、峰注釋、富集分析，下面將對各部分內(nèi)容進(jìn)行展開講解。
下機(jī)數(shù)據(jù)經(jīng)過過濾去除接頭含量過高或低質(zhì)量的reads，得到clean reads用于后續(xù)分析。常見的trim軟件有Trimmomatic、Skewer、fastp等。fastp是一款比較新的軟件，使用時可以用--adapter_sequence/--adapter_sequence_r2參數(shù)傳入接頭序列，也可以不填這兩個參數(shù)，軟件會自動識別接頭并進(jìn)行剪切。如：
fastp \
--in1 A1_1.fq.gz \ # read1原始fq文件
--out1 A1_clean_1.fq.gz \ # read1過濾后輸出的fq文件
--in2 A1_2.fq.gz ?\ # read2原始fq文件
--out2 A1_clean_2.fq.gz \ # read2過濾后輸出的fq文件
--cut_tail ?\ #從3’端向5’端滑窗，如果窗口內(nèi)堿基的平均質(zhì)量值小于設(shè)定閾值，則剪切
--cut_tail_window_size=1 \ #窗口大小
--cut_tail_mean_quality=30 \ #cut_tail參數(shù)對應(yīng)的平均質(zhì)量閾值
--average_qual=30 \ #如果一條read的堿基平均質(zhì)量值小于該值即會被舍棄
--length_required=20 ?\ #經(jīng)過剪切后的reads長度如果小于該值會被舍棄
fastp軟件的詳細(xì)使用方法可參考：。fastp軟件對于trim結(jié)果會生成網(wǎng)頁版的報告，可參考官網(wǎng)示例和，也可以用FastQC軟件對trim前后的數(shù)據(jù)質(zhì)量進(jìn)行評估，F(xiàn)astQC軟件會對單端的數(shù)據(jù)給出結(jié)果，如果是PE測序需要分別運(yùn)行兩次來評估read1和read2的數(shù)據(jù)質(zhì)量。
如：
fastqc A1_1.fq.gz
fastqc A1_2.fq.gz
FastQC會對reads從堿基質(zhì)量、接頭含量、N含量、高重復(fù)序列等多個方面對reads質(zhì)量進(jìn)行評估，生成詳細(xì)的網(wǎng)頁版報告，可參考官網(wǎng)示例：
經(jīng)過trim得到的reads可以使用BWA、bowtie2等軟件進(jìn)行比對。首先需要確定參考基因組fa文件，對fa文件建立索引。不同的軟件有各自建立索引的命令，BWA軟件可以參考如下方式建立索引：
bwa index genome.fa
建立好索引后即可開始比對，ATAC-seq推薦使用mem算法，輸出文件經(jīng)samtools排序輸出bam：
bwa mem genome.fa ?A1_clean_1.fq.gz A1_clean_2.fq.gz
| samtools sort -O bam -T A1 A1.bam
值得注意的是，在實(shí)驗(yàn)過程中質(zhì)體并不能完全去除，因此會有部分reads比對到質(zhì)體序列上，需要去除比對到質(zhì)體上的序列，去除質(zhì)體序列可以通過samtools提取，具體方法如下：首先將不含質(zhì)體的染色體名稱寫到一個chrlist文件中，一條染色體的名稱寫成一行，然后執(zhí)行如下命令即可得到去除質(zhì)體的bam
samtools view -b A1.bam $chrlist A1.del_MT_PT.bam
用于后續(xù)分析的reads需要時唯一比對且去重復(fù)的，bwa比對結(jié)果可以通過MAPQ值來提取唯一比對reads，可以用picard、sambamba等軟件去除dup，最終得到唯一比對且去重復(fù)的bam文件。
比對后得到的bam文件可以轉(zhuǎn)化為bigWig（bw）格式，通過可視化軟件進(jìn)行展示。deeptools軟件可以實(shí)現(xiàn)bw格式轉(zhuǎn)化和可視化展示。首先需要在linux環(huán)境中安裝deeptools軟件，可以用以下命令實(shí)現(xiàn)bam向bw格式的轉(zhuǎn)換：
bamCoverage -b A1.bam -o A1.bw
此外，可以使用deeptools軟件展示reads在特定區(qū)域的分布，如：
computeMatrix reference-point ??\ # reference-pioint表示計算一個參照點(diǎn)附近的reads分布，與之相對的是scale-regions，計算一個區(qū)域附近的reads分布
--referencePoint TSS ??\#以輸入的bed文件的起始位置作為參照點(diǎn)
-S ?A1.bw \ #可以是一個或多個bw文件
-R ?gene.bed \ #基因組位置文件
-b 3000 ??\ #計算邊界為參考點(diǎn)上游3000bp
-a 3000 ??\ #計算邊界為參考點(diǎn)下游3000bp，與-b合起來就是繪制參考點(diǎn)上下游3000bp以內(nèi)的reads分布
-o ?A1.matrix.mat.gz \ #輸出作圖數(shù)據(jù)名稱
#圖形繪制
plotHeatmap \
-m ?new_A1.matrix.mat.gz \ #上一步生成的作圖數(shù)據(jù)
-out A1.pdf \ # 輸出圖片名稱
繪圖結(jié)果展示：
MACS2能夠檢測DNA片斷的富集區(qū)域，是ATAC-seq數(shù)據(jù)call peak的主流軟件。峰檢出的原理如下：首先將所有的reads都向3'方向延伸插入片段長度，然后將基因組進(jìn)行滑窗，計算該窗口的dynamic λ，λ的計算公式為：λlocal = λBG（λBG是指背景區(qū)域上的reads數(shù)目），然后利用泊松分布模型的公式計算該窗口的顯著性P值，最后對每一個窗口的顯著性P值進(jìn)行FDR校正。默認(rèn)校正后的P值（即qvalue）小于或者等于0.05的區(qū)域?yàn)閜eak區(qū)域。需要現(xiàn)在linux環(huán)境中安裝macs2軟件，然后執(zhí)行以下命令：
macs2 callpeak \
-t A1.uni.dedup.bam \ #bam文件
-n A1 \ # 輸出文件前綴名
--shift -100 \ #extsize的一半乘以-1
--extsize 200 \ #一般是核小體大小
--call-summits #檢測峰頂信息
注：以上參數(shù)參考文獻(xiàn)（Jie Wang，et.al.2018.“ATAC-Seq analysis reveals a widespread decrease of chromatin accessibility in age-related macular degeneration.”Nature Communications）
ATAC分析得到的peak是染色質(zhì)上的開放區(qū)域，這些染色質(zhì)開放區(qū)域常常預(yù)示著轉(zhuǎn)錄因子的結(jié)合，因此對peak區(qū)域進(jìn)行motif分析很有意義。常見的motif分析軟件有homer和MEME。以homer軟件為例，首先在linux環(huán)境中安裝homer，然后用以下命令進(jìn)行motif分析：
findMotifsGenome.pl \
A1_peaks.bed \ #用于進(jìn)行motif分析的bed文件
genome.fa ?\ #參考基因組fa文件
A1 ?\ #輸出文件前綴
-size ?given \ #使用給定的bed區(qū)域位置進(jìn)行分析，如果填-size -100,50則是用給定bed中間位置的上游100bp到下游50bp的區(qū)域進(jìn)行分析
homer分析motif的原理及結(jié)果參見：
根據(jù)motif與已知轉(zhuǎn)錄因子的富集情況可以繪制氣泡圖，從而可以看到樣本與已知轉(zhuǎn)錄因子的富集顯著性。
差異peak代表著比較組合染色質(zhì)開放性有差異的位點(diǎn)，ChIP-seq和ATAC-seq都可以用DiffBind進(jìn)行差異分析。DiffBind通過可以通過bam文件和peak的bed文件計算出peak區(qū)域標(biāo)準(zhǔn)化的readcount，可以選擇edgeR、DESeq2等模型進(jìn)行差異分析。
在科研分析中我們往往需要將peak區(qū)域與基因聯(lián)系起來，也就是通過對peak進(jìn)行注釋找到peak相關(guān)基因。常見的peak注釋軟件有ChIPseeker、homer、PeakAnnotator等。以ChIPseeker為例，需要在R中安裝ChIPseeker包和GenomicFeatures包，然后就可以進(jìn)行分析了。
library(ChIPseeker)
library(GenomicFeatures)
txdb- makeTxDbFromGFF(‘gene.gtf’)#生成txdb對象，如果研究物種沒有已知的TxDb,可以用GenomicFeatures中的函數(shù)生成
peakfile -readPeakFile(‘A1_peaks.narrowPeak’)#導(dǎo)入需要注釋的peak文件
peakAnno - annotatePeak(peakfile,tssRegion=c(-2000, 2000), TxDb=txdb)
# 用peak文件和txdb進(jìn)行peak注釋，這里可以通過tssRegion定義TSS區(qū)域的區(qū)間
對于peak注釋的結(jié)果，也可以進(jìn)行可視化展示，如：
p - plotAnnoPie(peakAnno)
通過注釋得到的peak相關(guān)基因可以使用goseq、topGO等R包進(jìn)行GO富集分析，用kobas進(jìn)行kegg富集分析，也可以使用DAVID在線工具來完成富集分析。可以通過挑選感興趣的GO term或pathway進(jìn)一步篩選候選基因。
歡迎來到"bio生物信息"的世界
早期的研究普遍只做常染色體的全基因組關(guān)聯(lián)分析，很少做性染色體的。
主要原因是性染色體的遺傳模式比較復(fù)雜，存在X染色體失活，而且男女效應(yīng)值不大一樣。
其次，也不是所有的表型都是男女有差異的。
再然后，也沒有很好的工具計算性染色體的關(guān)聯(lián)分析。
隨著遺傳學(xué)的研究發(fā)展，現(xiàn)在有很多工具是允許計算性染色體的關(guān)聯(lián)分析。
下面簡單介紹一個常見的工具 SNPTEST
SNPTEST支持很多分析
比如，
對于linux系統(tǒng)而言，建議選擇動態(tài)鏈接版本（文件寫著dynamic）
wget
tar zxvf snptest_v2.5.4-beta3_CentOS6.6_x86_64_dynamic.tgz
輸入文件需要兩種類型。一種是表型文件，以 .sample 后綴，一種是基因型文件。
下圖是表型文件的格式
第一行是表型的title，第二行是對每一列的數(shù)據(jù)說明。
注意， 頭兩行是必須的 ，不然會報錯。
先講第一行的格式：
第一列和第二列是樣本的family ID 和個體ID。
第三列是missing，指的是樣本的缺失率，這一列可以通過plink的 --missing 參數(shù)獲得。
第四列到第七列都是協(xié)變量。(紅色框框)
第八列到第十一列都是表型。(藍(lán)色框框)
最后一列是性別。(綠色框框)
再講第二行的格式：
第二行的 0 0 0 D D C C P P B B D 又是什么呢
前三個 0 0 0 不需要修改，直接照著寫。
紅色框框 D D C C 指的是協(xié)變量的類型為離散型（D）和連續(xù)型（C）
藍(lán)色框框 P P B B 指的是表型的類型為連續(xù)型（P）和二分類（B）
綠色框框 D 指的是性別為離散型（D）
基因型文件支持三種格式。
第一種：GEN 或 gzipped GEN 格式，以.gen 或 .gen.gz結(jié)尾
第二種：BGEN格式，以.bgen結(jié)尾
第三種：plink格式，以.bed結(jié)尾
輸入如下命令：
./snptest \
-data ./example/cohort1_0X.bed ./example/cohort1.sample ./example/cohort2_0X.bed ./example/cohort2.sample \
-o ./example/ex.out \
-method newml \
-frequentist 1 \
-pheno bin1
解釋一下這些參數(shù)的意思。
-data 后面跟的是一個或多個隊列的基因型文件（.bed）和表型文件（.sample），這里列舉了兩個隊列。在實(shí)際的分析中，可以只分析一個，也可以同時分析多個隊列。
-o 指的是輸出的文件路徑（./example/）和文件名（ex.out）。
-method 指的是所用的方法。
-frequentist 指的是用的模型。模型可選加性模型、顯性模型、隱性模型、常規(guī)模型、雜合子模型。分別用1,2,3,4,5表示。 1=Additive, 2=Dominant, 3=Recessive, 4=General and 5=Heterozygote
-pheno 指的是所分析的表型列名。
報錯1：!! Error: (genfile::DuplicateIndividualError) A duplicate sample occurs on line 4 of the file
解決方法：這個報錯說明ID_1的字段是一樣的。需要將ID_1的每個樣本修改為獨(dú)一無二的字符?？梢耘cID_2保持一致。
報錯2：!! Error: the number of individuals (xxx) in the sample file differs from the number (yyy) in the genotypes file
解決方法：將基因型文件（.bed）的順序和數(shù)量與表型文件（.sample）的順序和數(shù)量保持一致
報錯3：二分類表型識別不了
解決方法：將二分類表型修改撐0,1編碼，SNPtest識別不了1,2
重慶分公司
重慶分公司
