拆分染色体linux命令 如何分离染色体

如何用vcftools从VCF文件中提取某条染色体信息

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

我们提供的服务有:成都做网站、网站建设、微信公众号开发、网站优化、网站认证、门头沟ssl等。为数千家企事业单位解决了网站和推广的问题。提供周到的售前咨询和贴心的售后服务,是有科学管理、有技术的门头沟网站制作公司

说明:

–gzvcf:处理压缩格式的vcf文件(可替换为–vcf)

–chr n:选择染色体n,例:–chr 1

–recode:重新编码为vcf文件,有过滤操作都要加上--recode

–recode-INFO-all:将输出的文件保存所有INFO信息

–stdout:标准输出,后接管道命令

–gzip -c:压缩

--max-missing

--max-missing的取值是0-1,为1时表示某个位点上所有的样本必须都有基因型,一个样本的基因型都不能缺。所以这个选项可以理解为:能分型的样本占总样本的比例至少为多少。

基本的思想就是利用数据流重定向,把原来输出到屏幕上的数据定向""到文件里

ATAC-seq专题---生信分析流程

ATAC-seq信息分析流程主要分为以下几个部分:数据质控、序列比对、峰检测、motif分析、峰注释、富集分析,下面将对各部分内容进行展开讲解。

下机数据经过过滤去除接头含量过高或低质量的reads,得到clean reads用于后续分析。常见的trim软件有Trimmomatic、Skewer、fastp等。fastp是一款比较新的软件,使用时可以用--adapter_sequence/--adapter_sequence_r2参数传入接头序列,也可以不填这两个参数,软件会自动识别接头并进行剪切。如:

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’端滑窗,如果窗口内碱基的平均质量值小于设定阈值,则剪切

--cut_tail_window_size=1 \ #窗口大小

--cut_tail_mean_quality=30 \ #cut_tail参数对应的平均质量阈值

--average_qual=30 \ #如果一条read的碱基平均质量值小于该值即会被舍弃

--length_required=20  \ #经过剪切后的reads长度如果小于该值会被舍弃

fastp软件的详细使用方法可参考:。fastp软件对于trim结果会生成网页版的报告,可参考官网示例和,也可以用FastQC软件对trim前后的数据质量进行评估,FastQC软件会对单端的数据给出结果,如果是PE测序需要分别运行两次来评估read1和read2的数据质量。

如:

fastqc A1_1.fq.gz

fastqc A1_2.fq.gz

FastQC会对reads从碱基质量、接头含量、N含量、高重复序列等多个方面对reads质量进行评估,生成详细的网页版报告,可参考官网示例:

经过trim得到的reads可以使用BWA、bowtie2等软件进行比对。首先需要确定参考基因组fa文件,对fa文件建立索引。不同的软件有各自建立索引的命令,BWA软件可以参考如下方式建立索引:

bwa index genome.fa

建立好索引后即可开始比对,ATAC-seq推荐使用mem算法,输出文件经samtools排序输出bam:

bwa mem genome.fa  A1_clean_1.fq.gz A1_clean_2.fq.gz

| samtools sort -O bam -T A1 A1.bam

值得注意的是,在实验过程中质体并不能完全去除,因此会有部分reads比对到质体序列上,需要去除比对到质体上的序列,去除质体序列可以通过samtools提取,具体方法如下:首先将不含质体的染色体名称写到一个chrlist文件中,一条染色体的名称写成一行,然后执行如下命令即可得到去除质体的bam

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

用于后续分析的reads需要时唯一比对且去重复的,bwa比对结果可以通过MAPQ值来提取唯一比对reads,可以用picard、sambamba等软件去除dup,最终得到唯一比对且去重复的bam文件。

比对后得到的bam文件可以转化为bigWig(bw)格式,通过可视化软件进行展示。deeptools软件可以实现bw格式转化和可视化展示。首先需要在linux环境中安装deeptools软件,可以用以下命令实现bam向bw格式的转换:

bamCoverage -b A1.bam -o A1.bw

此外,可以使用deeptools软件展示reads在特定区域的分布,如:

computeMatrix reference-point   \ # reference-pioint表示计算一个参照点附近的reads分布,与之相对的是scale-regions,计算一个区域附近的reads分布

--referencePoint TSS   \#以输入的bed文件的起始位置作为参照点

-S  A1.bw \ #可以是一个或多个bw文件

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

-b 3000   \ #计算边界为参考点上游3000bp

-a 3000   \ #计算边界为参考点下游3000bp,与-b合起来就是绘制参考点上下游3000bp以内的reads分布

-o  A1.matrix.mat.gz \ #输出作图数据名称

#图形绘制

plotHeatmap \

-m  new_A1.matrix.mat.gz \ #上一步生成的作图数据

-out A1.pdf \ # 输出图片名称

绘图结果展示:

MACS2能够检测DNA片断的富集区域,是ATAC-seq数据call peak的主流软件。峰检出的原理如下:首先将所有的reads都向3'方向延伸插入片段长度,然后将基因组进行滑窗,计算该窗口的dynamic λ,λ的计算公式为:λlocal = λBG(λBG是指背景区域上的reads数目),然后利用泊松分布模型的公式计算该窗口的显著性P值,最后对每一个窗口的显著性P值进行FDR校正。默认校正后的P值(即qvalue)小于或者等于0.05的区域为peak区域。需要现在linux环境中安装macs2软件,然后执行以下命令:

macs2 callpeak \

-t A1.uni.dedup.bam \ #bam文件

-n A1 \ # 输出文件前缀名

--shift -100 \ #extsize的一半乘以-1

--extsize 200 \ #一般是核小体大小

--call-summits #检测峰顶信息

注:以上参数参考文献(Jie Wang,et.al.2018.“ATAC-Seq analysis reveals a widespread decrease of chromatin accessibility in age-related macular degeneration.”Nature Communications)

ATAC分析得到的peak是染色质上的开放区域,这些染色质开放区域常常预示着转录因子的结合,因此对peak区域进行motif分析很有意义。常见的motif分析软件有homer和MEME。以homer软件为例,首先在linux环境中安装homer,然后用以下命令进行motif分析:

findMotifsGenome.pl \

A1_peaks.bed \ #用于进行motif分析的bed文件

genome.fa  \ #参考基因组fa文件

A1  \ #输出文件前缀

-size  given \ #使用给定的bed区域位置进行分析,如果填-size -100,50则是用给定bed中间位置的上游100bp到下游50bp的区域进行分析

homer分析motif的原理及结果参见:

根据motif与已知转录因子的富集情况可以绘制气泡图,从而可以看到样本与已知转录因子的富集显著性。

差异peak代表着比较组合染色质开放性有差异的位点,ChIP-seq和ATAC-seq都可以用DiffBind进行差异分析。DiffBind通过可以通过bam文件和peak的bed文件计算出peak区域标准化的readcount,可以选择edgeR、DESeq2等模型进行差异分析。

在科研分析中我们往往需要将peak区域与基因联系起来,也就是通过对peak进行注释找到peak相关基因。常见的peak注释软件有ChIPseeker、homer、PeakAnnotator等。以ChIPseeker为例,需要在R中安装ChIPseeker包和GenomicFeatures包,然后就可以进行分析了。

library(ChIPseeker)

library(GenomicFeatures)

txdb- makeTxDbFromGFF(‘gene.gtf’)#生成txdb对象,如果研究物种没有已知的TxDb,可以用GenomicFeatures中的函数生成

peakfile -readPeakFile(‘A1_peaks.narrowPeak’)#导入需要注释的peak文件

peakAnno - annotatePeak(peakfile,tssRegion=c(-2000, 2000), TxDb=txdb)

# 用peak文件和txdb进行peak注释,这里可以通过tssRegion定义TSS区域的区间

对于peak注释的结果,也可以进行可视化展示,如:

p - plotAnnoPie(peakAnno)

通过注释得到的peak相关基因可以使用goseq、topGO等R包进行GO富集分析,用kobas进行kegg富集分析,也可以使用DAVID在线工具来完成富集分析。可以通过挑选感兴趣的GO term或pathway进一步筛选候选基因。

对性染色体进行关联分析

欢迎来到"bio生物信息"的世界

早期的研究普遍只做常染色体的全基因组关联分析,很少做性染色体的。

主要原因是性染色体的遗传模式比较复杂,存在X染色体失活,而且男女效应值不大一样。

其次,也不是所有的表型都是男女有差异的。

再然后,也没有很好的工具计算性染色体的关联分析。

随着遗传学的研究发展,现在有很多工具是允许计算性染色体的关联分析。

下面简单介绍一个常见的工具 SNPTEST

SNPTEST支持很多分析

比如,

对于linux系统而言,建议选择动态链接版本(文件写着dynamic)

wget

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

输入文件需要两种类型。一种是表型文件,以 .sample 后缀,一种是基因型文件。

下图是表型文件的格式

第一行是表型的title,第二行是对每一列的数据说明。

注意, 头两行是必须的 ,不然会报错。

先讲第一行的格式:

第一列和第二列是样本的family ID 和个体ID。

第三列是missing,指的是样本的缺失率,这一列可以通过plink的 --missing 参数获得。

第四列到第七列都是协变量。(红色框框)

第八列到第十一列都是表型。(蓝色框框)

最后一列是性别。(绿色框框)

再讲第二行的格式:

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

前三个 0 0 0 不需要修改,直接照着写。

红色框框 D D C C 指的是协变量的类型为离散型(D)和连续型(C)

蓝色框框 P P B B 指的是表型的类型为连续型(P)和二分类(B)

绿色框框 D 指的是性别为离散型(D)

基因型文件支持三种格式。

第一种:GEN 或 gzipped GEN 格式,以.gen 或 .gen.gz结尾

第二种:BGEN格式,以.bgen结尾

第三种:plink格式,以.bed结尾

输入如下命令:

./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

解释一下这些参数的意思。

-data 后面跟的是一个或多个队列的基因型文件(.bed)和表型文件(.sample),这里列举了两个队列。在实际的分析中,可以只分析一个,也可以同时分析多个队列。

-o 指的是输出的文件路径(./example/)和文件名(ex.out)。

-method 指的是所用的方法。

-frequentist 指的是用的模型。模型可选加性模型、显性模型、隐性模型、常规模型、杂合子模型。分别用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的每个样本修改为独一无二的字符。可以与ID_2保持一致。

报错2:!! Error: the number of individuals (xxx) in the sample file differs from the number (yyy) in the genotypes file

解决方法:将基因型文件(.bed)的顺序和数量与表型文件(.sample)的顺序和数量保持一致

报错3:二分类表型识别不了

解决方法:将二分类表型修改撑0,1编码,SNPtest识别不了1,2

当前名称:拆分染色体linux命令 如何分离染色体
转载源于:https://www.cdcxhl.com/article40/doicoho.html

成都网站建设公司_创新互联,为您提供网站设计公司关键词优化搜索引擎优化网站改版品牌网站制作App开发

广告

声明:本网站发布的内容(图片、视频和文字)以用户投稿、用户转载内容为主,如果涉及侵权请尽快告知,我们将会在第一时间删除。文章观点不代表本网站立场,如需处理请联系客服。电话:028-86922220;邮箱:631063699@qq.com。内容未经允许不得转载,或转载时需注明来源: 创新互联

成都定制网站网页设计