拆分染色體linux命令 如何分離染色體

如何用vcftools從VCF文件中提取某條染色體信息

vcftools --gzvcf input.vcf --chr n --recode – recode-INFO-all --stdout | gzip -c output.vcf.gz

我們提供的服務(wù)有：成都做網(wǎng)站、網(wǎng)站建設(shè)、微信公眾號(hào)開發(fā)、網(wǎng)站優(yōu)化、網(wǎng)站認(rèn)證、門頭溝ssl等。為數(shù)千家企事業(yè)單位解決了網(wǎng)站和推廣的問(wèn)題。提供周到的售前咨詢和貼心的售后服務(wù)，是有科學(xué)管理、有技術(shù)的門頭溝網(wǎng)站制作公司

說(shuō)明：

–gzvcf：處理壓縮格式的vcf文件（可替換為–vcf）

–chr n：選擇染色體n，例：–chr 1

–recode：重新編碼為vcf文件,有過(guò)濾操作都要加上--recode

–recode-INFO-all：將輸出的文件保存所有INFO信息

–stdout：標(biāo)準(zhǔn)輸出，后接管道命令

–gzip -c：壓縮

--max-missing

--max-missing的取值是0-1，為1時(shí)表示某個(gè)位點(diǎn)上所有的樣本必須都有基因型，一個(gè)樣本的基因型都不能缺。所以這個(gè)選項(xiàng)可以理解為：能分型的樣本占總樣本的比例至少為多少。

基本的思想就是利用數(shù)據(jù)流重定向，把原來(lái)輸出到屏幕上的數(shù)據(jù)定向""到文件里

ATAC-seq專題---生信分析流程

ATAC-seq信息分析流程主要分為以下幾個(gè)部分：數(shù)據(jù)質(zhì)控、序列比對(duì)、峰檢測(cè)、motif分析、峰注釋、富集分析，下面將對(duì)各部分內(nèi)容進(jìn)行展開講解。

下機(jī)數(shù)據(jù)經(jīng)過(guò)過(guò)濾去除接頭含量過(guò)高或低質(zhì)量的reads，得到clean reads用于后續(xù)分析。常見的trim軟件有Trimmomatic、Skewer、fastp等。fastp是一款比較新的軟件，使用時(shí)可以用--adapter_sequence/--adapter_sequence_r2參數(shù)傳入接頭序列，也可以不填這兩個(gè)參數(shù)，軟件會(huì)自動(dòng)識(shí)別接頭并進(jìn)行剪切。如：

fastp \

--in1 A1_1.fq.gz \ # read1原始fq文件

--out1 A1_clean_1.fq.gz \ # read1過(guò)濾后輸出的fq文件

--in2 A1_2.fq.gz ?\ # read2原始fq文件

--out2 A1_clean_2.fq.gz \ # read2過(guò)濾后輸出的fq文件

--cut_tail ?\ #從3’端向5’端滑窗，如果窗口內(nèi)堿基的平均質(zhì)量值小于設(shè)定閾值，則剪切

--cut_tail_window_size=1 \ #窗口大小

--cut_tail_mean_quality=30 \ #cut_tail參數(shù)對(duì)應(yīng)的平均質(zhì)量閾值

--average_qual=30 \ #如果一條read的堿基平均質(zhì)量值小于該值即會(huì)被舍棄

--length_required=20 ?\ #經(jīng)過(guò)剪切后的reads長(zhǎng)度如果小于該值會(huì)被舍棄

fastp軟件的詳細(xì)使用方法可參考：。fastp軟件對(duì)于trim結(jié)果會(huì)生成網(wǎng)頁(yè)版的報(bào)告，可參考官網(wǎng)示例和，也可以用FastQC軟件對(duì)trim前后的數(shù)據(jù)質(zhì)量進(jìn)行評(píng)估，F(xiàn)astQC軟件會(huì)對(duì)單端的數(shù)據(jù)給出結(jié)果，如果是PE測(cè)序需要分別運(yùn)行兩次來(lái)評(píng)估read1和read2的數(shù)據(jù)質(zhì)量。

如：

fastqc A1_1.fq.gz

fastqc A1_2.fq.gz

FastQC會(huì)對(duì)reads從堿基質(zhì)量、接頭含量、N含量、高重復(fù)序列等多個(gè)方面對(duì)reads質(zhì)量進(jìn)行評(píng)估，生成詳細(xì)的網(wǎng)頁(yè)版報(bào)告，可參考官網(wǎng)示例：

經(jīng)過(guò)trim得到的reads可以使用BWA、bowtie2等軟件進(jìn)行比對(duì)。首先需要確定參考基因組fa文件，對(duì)fa文件建立索引。不同的軟件有各自建立索引的命令，BWA軟件可以參考如下方式建立索引：

bwa index genome.fa

建立好索引后即可開始比對(duì)，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)過(guò)程中質(zhì)體并不能完全去除，因此會(huì)有部分reads比對(duì)到質(zhì)體序列上，需要去除比對(duì)到質(zhì)體上的序列，去除質(zhì)體序列可以通過(guò)samtools提取，具體方法如下：首先將不含質(zhì)體的染色體名稱寫到一個(gè)chrlist文件中，一條染色體的名稱寫成一行，然后執(zhí)行如下命令即可得到去除質(zhì)體的bam

samtools view -b A1.bam $chrlist A1.del_MT_PT.bam

用于后續(xù)分析的reads需要時(shí)唯一比對(duì)且去重復(fù)的，bwa比對(duì)結(jié)果可以通過(guò)MAPQ值來(lái)提取唯一比對(duì)reads，可以用picard、sambamba等軟件去除dup，最終得到唯一比對(duì)且去重復(fù)的bam文件。

比對(duì)后得到的bam文件可以轉(zhuǎn)化為bigWig（bw）格式，通過(guò)可視化軟件進(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表示計(jì)算一個(gè)參照點(diǎn)附近的reads分布，與之相對(duì)的是scale-regions，計(jì)算一個(gè)區(qū)域附近的reads分布

--referencePoint TSS ??\#以輸入的bed文件的起始位置作為參照點(diǎn)

-S ?A1.bw \ #可以是一個(gè)或多個(gè)bw文件

-R ?gene.bed \ #基因組位置文件

-b 3000 ??\ #計(jì)算邊界為參考點(diǎn)上游3000bp

-a 3000 ??\ #計(jì)算邊界為參考點(diǎn)下游3000bp，與-b合起來(lái)就是繪制參考點(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能夠檢測(cè)DNA片斷的富集區(qū)域，是ATAC-seq數(shù)據(jù)call peak的主流軟件。峰檢出的原理如下：首先將所有的reads都向3'方向延伸插入片段長(zhǎng)度，然后將基因組進(jìn)行滑窗，計(jì)算該窗口的dynamic λ，λ的計(jì)算公式為：λlocal = λBG（λBG是指背景區(qū)域上的reads數(shù)目），然后利用泊松分布模型的公式計(jì)算該窗口的顯著性P值，最后對(duì)每一個(gè)窗口的顯著性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 #檢測(cè)峰頂信息

注：以上參數(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é)合，因此對(duì)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通過(guò)可以通過(guò)bam文件和peak的bed文件計(jì)算出peak區(qū)域標(biāo)準(zhǔn)化的readcount，可以選擇edgeR、DESeq2等模型進(jìn)行差異分析。

在科研分析中我們往往需要將peak區(qū)域與基因聯(lián)系起來(lái)，也就是通過(guò)對(duì)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對(duì)象，如果研究物種沒有已知的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注釋，這里可以通過(guò)tssRegion定義TSS區(qū)域的區(qū)間

對(duì)于peak注釋的結(jié)果，也可以進(jìn)行可視化展示，如：

p - plotAnnoPie(peakAnno)

通過(guò)注釋得到的peak相關(guān)基因可以使用goseq、topGO等R包進(jìn)行GO富集分析，用kobas進(jìn)行kegg富集分析，也可以使用DAVID在線工具來(lái)完成富集分析。可以通過(guò)挑選感興趣的GO term或pathway進(jìn)一步篩選候選基因。

對(duì)性染色體進(jìn)行關(guān)聯(lián)分析

歡迎來(lái)到"bio生物信息"的世界

早期的研究普遍只做常染色體的全基因組關(guān)聯(lián)分析，很少做性染色體的。

主要原因是性染色體的遺傳模式比較復(fù)雜，存在X染色體失活，而且男女效應(yīng)值不大一樣。

其次，也不是所有的表型都是男女有差異的。

再然后，也沒有很好的工具計(jì)算性染色體的關(guān)聯(lián)分析。

隨著遺傳學(xué)的研究發(fā)展，現(xiàn)在有很多工具是允許計(jì)算性染色體的關(guān)聯(lián)分析。

下面簡(jiǎn)單介紹一個(gè)常見的工具 SNPTEST

SNPTEST支持很多分析

比如，

對(duì)于linux系統(tǒng)而言，建議選擇動(dòng)態(tài)鏈接版本（文件寫著dynamic）

wget

tar zxvf snptest_v2.5.4-beta3_CentOS6.6_x86_64_dynamic.tgz

輸入文件需要兩種類型。一種是表型文件，以 .sample 后綴，一種是基因型文件。

下圖是表型文件的格式

第一行是表型的title，第二行是對(duì)每一列的數(shù)據(jù)說(shuō)明。

注意， 頭兩行是必須的 ，不然會(huì)報(bào)錯(cuò)。

先講第一行的格式：

第一列和第二列是樣本的family ID 和個(gè)體ID。

第三列是missing，指的是樣本的缺失率，這一列可以通過(guò)plink的 --missing 參數(shù)獲得。

第四列到第七列都是協(xié)變量。(紅色框框)

第八列到第十一列都是表型。(藍(lán)色框框)

最后一列是性別。(綠色框框)

再講第二行的格式：

第二行的 0 0 0 D D C C P P B B D 又是什么呢

前三個(gè) 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 后面跟的是一個(gè)或多個(gè)隊(duì)列的基因型文件（.bed）和表型文件（.sample），這里列舉了兩個(gè)隊(duì)列。在實(shí)際的分析中，可以只分析一個(gè)，也可以同時(shí)分析多個(gè)隊(duì)列。

-o 指的是輸出的文件路徑（./example/）和文件名（ex.out）。

-method 指的是所用的方法。

-frequentist 指的是用的模型。模型可選加性模型、顯性模型、隱性模型、常規(guī)模型、雜合子模型。分別用1,2,3,4,5表示。 1=Additive, 2=Dominant, 3=Recessive, 4=General and 5=Heterozygote

-pheno 指的是所分析的表型列名。

報(bào)錯(cuò)1：!! Error: (genfile::DuplicateIndividualError) A duplicate sample occurs on line 4 of the file

解決方法：這個(gè)報(bào)錯(cuò)說(shuō)明ID_1的字段是一樣的。需要將ID_1的每個(gè)樣本修改為獨(dú)一無(wú)二的字符。可以與ID_2保持一致。

報(bào)錯(cuò)2：!! Error: the number of individuals (xxx) in the sample file differs from the number (yyy) in the genotypes file

解決方法：將基因型文件（.bed）的順序和數(shù)量與表型文件（.sample）的順序和數(shù)量保持一致

報(bào)錯(cuò)3：二分類表型識(shí)別不了

解決方法：將二分類表型修改撐0,1編碼，SNPtest識(shí)別不了1,2

新聞名稱：拆分染色體linux命令 如何分離染色體
文章起源：http://www.jiaotiyi.com/article/doicoho.html