转录组分析流程

一、读文章获取下载数据

1、读文章

一般我都从NCBI上面下载文章，找到数据号

2、下载数据

进入NCBI的GEO数据库，输入数据号，通过FTP下载数据，并使用fastq-dump命令将SRR格式转换为FASTQ格式

for ((i=56;i<=62;i++));

do

fastq-dump --gzip --split-files SRR35899$i.sra

done 

二、质控

（1）、fastqc -o ../fastqc *.fastq.gz

（2）、数据过滤 trimmomatic;fastx_toolkit

三、下载参考基因组及注释序列（无参基因组需要拼接转录本TRINITY）

一般ensembl数据库下载FATSTA及GTF文件

ftp://ftp.ensembl.org/pub/

四、序列比对

使用HISAT2

（1）、构建索引序列

~/Python-2.7.14/python ~/software/hisat2-2.1.0/extract_exons.py Homo_sapiens.GRCh38.83.chr.gtf>genome.exon

~/Python-2.7.14/python ~/software/hisat2-2.1.0/
extract_splice_sites.py Homo_sapiens.GRCh38.83.chr.gtf>genome.ss

hisat2-build -p 32 Mus_musculus.GRCm38.dna.primary_assembly.fa --ss genome.ss --exon genome.exon genome_tran

（2）、序列比对

nohup hisat2 -p 16 --dta -x ~/20171014/refs/grch38/genome -1 ~/20171014/raw_data/SRR3589956_1.fastq.gz -2 ~/20171014/raw_data/SRR3589956_2.fastq.gz -S SRR3589956.sam

五、格式转换

samtools view -S    -b

samtools sort  -o

samtools index

六、计数

使用htseq进行计数

/stor9000/apps/users/NWSUAF/2016050428/Python-2.7.14/python  ~/Python-2.7.14/bin/htseq-count ~/20171014/hisat2/SRR3589958_count.sam ~/20171014/human/Homo_sapiens.GRCh38.90.chr.gtf

七、差异基因分析

（1）读取数据文件

control1<-read.table("SRR3589959.count",sep="\t",col.names=c("gene_id",control1))

control12<-read.table("SRR3589960.count",sep="\t",col.names=c("gene_id",control2))

rep1<-read.table("SRR3589961.count",sep="\t",col.names=c("gene_id",akap951))

rep2<-read.table("SRR3589962.count",sep="\t",col.names=c("gene_id",akap952))

raw_count<-merge(merge(control1,control2,by="gene_id),merge(rep1,rep2,by="gene_id"))

ensemble<-gsub("\\.\\*d","",raw_count$gene_id)

row.names(raw_count)=ensemble

raw_count_filter<-raw_count[,-1]

(2)构建DDS对象

condition<-factor(c(rep("control",2),rep("apak95",2)),levels=c("control","apak95"))

cound_data<-raw_count_filter[,1:4]

cc_data<-data.frame(row.name=col.names(raw_count_filter),condition)

dds<-DEseqDataSetFromMatrix(count_data,col_data,design=~condition)

(3)Deseq标准化dds

dds2<-Deseq(dds)

resultsname(dds2)

res<-result(dds2)

summary(res)

(4)提取差异分析结果

res<-res[order(res$padj),])

diff_gene_deseq2<-subset(res,padj<0.05 & (log2foldchange>1 | log2foldchange<-1))

resdata<merge(as.data.frame(res),as.data.frame(counts(dds2,normalize=TRUE)),by="row.names",sort=FALSE)

write.csv(as.data.frame(resdata),"差异分析结果")

八、GO、PATHWAY分析

enrichGO<-enrichGO(gene, OrgDb, keyType = "ENTREZID", ont = "MF",

  pvalueCutoff = 0.05, pAdjustMethod = "BH", universe, qvalueCutoff = 0.2,

  minGSSize = 10, maxGSSize = 500, readable = FALSE, pool = FALSE)

enrichPATHWAY<-enrichKEGG(gene, organism = "hsa", keyType = "kegg", pvalueCutoff = 0.05,

  pAdjustMethod = "BH", universe, minGSSize = 10, maxGSSize = 500,

  qvalueCutoff = 0.2, use_internal_data = FALSE)