基因家族分析套路（四）

ndsky 2017-06-05

展开全文

广告：大家有基因组家族分析的需求，麻烦留言哦，小编会第一时间回复的。

今天是基因家族分析类文章最后的一部分，也是一个文章亮点所在的部分，小伙伴们仔细阅读学习吧！

一、转录组及芯片原始数据下载网站

1、 GEO datesets/profile(http://www.ncbi.nlm./gds ).。

用法见下图。GEO数据ID命名规则：GPL->GSE->GSM.

GPL: platform

GSE: multiple series.

GSM: multiple samples.

GDS ≈ GSE. Thedifference concentrated on the data labeled GDS can be analyzed for one geneonline. It is simple and easily.

The data in the sameGPL can be used to compare inexperiment.

下面是在线分析转录组数据的用法：

2、EBI ArrayExpress(http://www./arrayexpress/)

该数据库下载数据用法如下：

3、PLEXdb(http://www./).

该数据库下载数据用法如下，注意用户名和密码！

4、SRA db(http://www.ncbi.nlm./sra/)

5、DRA db（http://trace.ddbj./DRASearch/）

二、数据处理

拿到原始数据，要进行处理，才能进行后续数据分析。

1、芯片数据。原始数据格式“.cel”格式。以AffyMicroarray数据处理为例讲述主要的命令如下：

> library(affy);

>library(makecdfenv);

>library……

> barleyGenome = make.cdf.env(“barleyGenome.cdf')

>mydata <- ReadAffy() ##choose “.cel “ file analyzed.

>eset <- rma(mydata);

>write.exprs(eset,file='mydata.txt')

>design <- model.matrix(~-1 factor(c(1,1,2,2,3,3))) # Createsappropriate design matrix.

>colnames(design) <-c('group1', 'group2', 'group3') # Assigns column names.

>fit <- lmFit(eset, design) # Fits a linear model for each gene based onthe given series of arrays.

>contrast.matrix <- makeContrasts(group2-group1,group3-group2, group3-group1, levels=design) # Creates appropriate contrast matrix toperform all pairwise comparisons.

>fit2 <- contrasts.fit(fit, contrast.matrix)# Computes estimatedcoefficients and standard errors for a given set of contrasts.

>fit2 <- eBayes(fit2) # Computes moderated t-statistics and log-oddsof differential expression by empirical Bayes

>topTable(fit2, coef=1,adjust='fdr', sort.by='B', number=10) # Generates list of top 10 ('number=10')differentially expressed genes sorted by B-values ('sort.by=B') for firstcomparison group.

>write.table(topTable(fit2, coef=1,adjust='fdr', sort.by='B', number=500),file='limma_complete.xls', row.names=F, sep='\t') # Exports complete limma statistics table forfirst comparison group.

>results <- decideTests(fit2,p.value=0.05); vennDiagram(results)

2、转录组数据处理。原始数据格式为sra或fastq格式。Sra可以转换为fastq然后运用下面的命令进行处理。

1）获得cleandata；

fastx_clipper :clip adapter.

fastq_quality_filter: base quality control.

fastq_quality_trimmer: trim 5’ low quality bases.

2）计算RPKM.

bowtie2-buildpath/db.seq path/db

tophat db read.fastq

bam_filter path/accepted_hits.bam

samtools view -h -o output-uniq.sam output_uniq.bam

excel for calculation(low frequencyreads ≤5 were omitted ).