一、读文章获取下载数据
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)
