- 數據格式如下
gene_id Sham-1 Sham-2 Sham-3 Sham-4 Sham-5 Rep-1h-1 Rep-1h-2 Rep-1h-3 Rep-1h-4 Rep-1h-5 Rep-3h-1 Rep-3h-2 Rep-3h-3 Rep-3h-4 Rep-3h-5 Rep-6h-1 Rep-6h-2 Rep-6h-3 Rep-6h-4 Rep-6h-5 Rep-12h-1 Rep-12h-2 Rep-12h-3 Rep-12h-4 Rep-12h-5 Rep-24h-1 Rep-24h-2 Rep-24h-3 Rep-24h-4 Rep-24h-5
Traf1 0.204531 0.358811 0.24649 0.270231 0.169885 0.588808 0.526418 0.87557 0.403861 1.846186 1.890555 1.041459 0.881003 1.036722 1.016675 2.09069 3.09231 1.180259 2.610673 0.888904 2.936677 1.829962 1.857749 2.93743 2.424644 1.3602 2.654057 2.121849 2.309851 2.052516
- 代碼如下：
a=read.table(file="Traf1.fpkm.lst")
c=t.data.frame(a)
gene_id=c(c[,1])
gene=as.vector(gene_id)
id=gene[2:length(gene)]
va=as.vector(c[,2])[2:length(as.vector(c[,2]))]
fpkm=round(as.numeric(va),3)
dd=data.frame(n,fpkm)
ggplot(dd,aes(x=n,y=fpkm))+
geom_bar(stat = "identity",fill ="#09BFFE",colour="grey")+
xlab("Samples")+ylab("FPKM of Traf1 in different sample")+
theme_classic()+
theme(axis.text.y=element_blank(),axis.text.x=element_text(size=8,color="black",angle = 30,face="bold",hjust = 1))+
geom_text(aes(label=dd$fpkm), position=position_stack(vjust=1.05),angle=30)+
scale_y_continuous(breaks = seq(0,3.5,0.5),limits = c(0,3.5),labels = seq(0,3.5,0.5))
使用fimo對Traf1基因進行motif的尋找
- 下載motif數據庫
wget http://meme-suite.org/meme-software/Databases/motifs/motif_databases.12.18.tgz
tar zxvf motif_databases.12.18.tgz
代碼如下
ln -s /media/sda/user/chengxu/project/cardiac-IR/lncRNA_sequencing/expression_data/mRNA_symbol_FPKM_matrix.txt ./
les mRNA_symbol_FPKM_matrix.txt|perl -F"\t" -lane 'print join("\t",@F[0..5])."\t".join("\t",@F[43..67])' |head -1 >Traf1.fpkm.lst
les mRNA_symbol_FPKM_matrix.txt|perl -F"\t" -lane 'print join("\t",@F[0..5])."\t".join("\t",@F[43..67])' |grep Traf1 >>Traf1.fpkm.lst
grep Traf1 /media/sda/database/Ensembl/release-89/mus_musculus/gtf/gene.bed.txt |perl -lane '$s=$F[5]-1;$e=$F[6]-1;print qq{$F[3]\t$s\t$e\t$F[7]}' >Traf1.bed
bedtools getfasta -fi /media/sda/database/Ensembl/release-89/mus_musculus/fasta/mm10.fasta -bed Traf1.bed -fo Traf1.gene.fasta
fimo -oc Traf1 JASPAR2018_CORE_vertebrates_non-redundant.meme Traf1.gene.fasta
- 需要注意的是bed文件時從0開始算位置的，而gtf等注釋文件時從1開始算的，所以gtf轉換成bed的時候要減少一個鹼基的位置
- 數據格式
aaa sss fpkm
Sham Sham-1 0.204531
Sham Sham-2 0.358811
Sham Sham-3 0.24649
Sham Sham-4 0.270231
Sham Sham-5 0.169885
Rep-1h Rep-1h-1 0.588808
Rep-1h Rep-1h-2 0.526418
Rep-1h Rep-1h-3 0.87557
Rep-1h Rep-1h-4 0.403861
Rep-1h Rep-1h-5 1.846186
Rep-3h Rep-3h-1 1.890555
Rep-3h Rep-3h-2 1.041459
Rep-3h Rep-3h-3 0.881003
Rep-3h Rep-3h-4 1.036722
Rep-3h Rep-3h-5 1.016675
Rep-6h Rep-6h-1 2.09069
Rep-6h Rep-6h-2 3.09231
Rep-6h Rep-6h-3 1.180259
Rep-6h Rep-6h-4 2.610673
Rep-6h Rep-6h-5 0.888904
Rep-12h Rep-12h-1 2.936677
Rep-12h Rep-12h-2 1.829962
Rep-12h Rep-12h-3 1.857749
Rep-12h Rep-12h-4 2.93743
Rep-12h Rep-12h-5 2.424644
Rep-24h Rep-24h-1 1.3602
Rep-24h Rep-24h-2 2.654057
Rep-24h Rep-24h-3 2.121849
Rep-24h Rep-24h-4 2.309851
Rep-24h Rep-24h-5 2.052516
- 畫圖代碼
dd=read.table(file = "Traf1.changeformat",header = T)
#修改箱線圖的順序
dd$aaa=factor(dd$aaa,levels=c("Sham","Rep-3h","Rep-6h","Rep-12h","Rep-24h"))
ggplot(dd, aes(x=aaa, y=fpkm, fill=aaa)) +
geom_boxplot() +
scale_fill_manual(values=c("#E41A1C","#377EB8","#4DAF4A","#984EA3","#FF7F00","#FFFF33"), name="Var2") +
labs(x="", y="FPKM of Traf1") + #去掉x軸的Var2
guides(fill = guide_legend(title="Group", keywidth=2, keyheight=2)) +
theme_bw() +
theme(panel.grid = element_blank()) +
theme(axis.title.y = element_text(size=18)) +
theme(axis.text.y = element_text(size=16, hjust=1)) +
theme(axis.text.x = element_text(colour="grey20", size=96/len, angle=30, hjust=1)) +
theme(legend.title = element_text(size=15)) +
theme(legend.text = element_text(size=15))
