linux生信分析?生信分析公司
今天给各位分享linux生信分析的知识,其中也会对生信分析公司进行解释,如果能碰巧解决你现在面临的问题,别忘了关注本站,现在开始吧!
用于生物信息分析该如何安装ubuntu系统
1.生信软件系统的选择——Linux(ubuntu)
对于生信分析人员来说,日常工作,软件运行,跑流程,均在linux下操作。当然,也有基于云端的生信分析平台,如免费的Galaxy,或者某些公司的一站式云平台。
比较初学者学生物信息还是使用开源软件、学原理、一步一步运行才有意思。这路子,一定要适应Linux的命令行界面。
选择windows还是linux?一定是linux,windows太多的生物软件不兼容了。
选择linux的哪个版本?推荐桌面版的Ubuntu——稳定,美观,适合初学者之称;次之,Centos——免费、稳定的服务器linux版本之称。
用那种方式安装linux好?推荐虚拟机安装。不太建议双系统,云端这种。因为,对于初学者在系统中,需要反复折腾,测试,搞垮系统是常事。
选择开源的VMbox还是商业版VMware?两者都可以,但各有缺点。VMbox更新比较快,经常更新后,可能会出现报错,系统无法打开的现象,较低版本的反而比较稳定,如果用好了,不建议经常更新。还有一点是,VMbox在鼠标控制上,没有VMware流畅。VMware十分稳定,流程好用。最新版一般要收费。可以选择比最新版版本稍低的,上网搜注册码,免费使用。还是那样,用好了,不要经常更新。某些生信软件会提供VMbox的镜像,如qiime。
VMbox的镜像能不能转到VMware上使用?,答案是可以的,使用VMbox的镜像导出功能,然后使用VMware进行导入,保持两者格式相同。
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进一步筛选候选基因。
生信分析怎么学
学习生信分析需要具备一定的计算机、生物学和统计学知识,建议按以下步骤学习:
1.建立基础知识:先学习生物学、计算机科学和统计学的基础知识,掌握常用的生物学术语和基本的编程概念。可以参考一些经典教材如《生物信息学导论》、《R语言实战》等。
2.学习常用工具和软件:学习生物信息学分析中常用的工具和软件,例如NCBI、BLAST、UCSC等数据库和软件,学习Linux操作系统和常用命令,掌握编程语言如Perl、Python、R等的使用。
3.参加课程或培训:参加一些线上或线下的课程或培训,例如Coursera上的生物信息学课程、培训班、讲座等,了解生物信息学分析的流程和方法,掌握实践技能。
4.实践和练习:通过实际项目的实践,积累经验和技能。可以通过参加竞赛、学术项目或者开源社区的项目来进行实践。
5.学习交流:通过参加学术会议、讨论组、社区等渠道,与其他从业人员交流和分享经验,了解最新的技术发展和应用实践。
总之,星科SCIER认为学习生信分析需要综合多个学科知识,需要不断实践和练习,才能熟练掌握相关技能。