RNA-seq生命科学-简书专题RNA Seq流程

RNA-seq分析文章结果重现

2018-07-26  本文已影响42人  Dawn_WangTP

文章：Temporal dynamics of gene expression and histone marks at the Arabidopsis shoot meristem during flowering.

一、查看并下载数据

tr 替换\t 为\n ； nl为对每行加行号
head -n1 E-MTAB-5130.sdrf.txt | tr '\t' '\n' | nl | grep URI
image 

tail -n +2 E-MTAB-5130.sdrf.txt | cut -f 33 | xargs -i wget {}

二、fastqc和MultiQC质控

/usr/local/bin/fastqc -o ./ --nogroup /sas/supercloud-kong/wangtianpeng/data/ath_flowering/1_raw_data/ERR1698201_2.fastq.gz 

### multiqc质控
multiqc *fastqc.zip --pdf

三、Trimmo修建

for i in `seq 194 209`; do  i=ERR1698$i; nohup trimmomatic PE ../1_raw_data/$i\_1.fastq.gz ../1_raw_data/$i\_2.fastq.gz $i\_P1.fastq.gz $i\_U1.fastq.gz $i\_P2.fastq.gz $i\_U2.fastq.gz ILLUMINACLIP:/home/wangtianpeng/anaconda3/share/trimmomatic-0.36-5/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:5 TRAILING:5 SLIDINGWINDOW:4:5 MINLEN:20 \&; done

四、SALMON定量

salmon index -t athal.fa.gz -i athal_index

for i in `seq 194 209`;do nohup salmon quant -i /sas/supercloud-kong/wangtianpeng/Database/plant/a_th/ath_index -l A -1 /sas/supercloud-kong/wangtianpeng/data/ath_flowering/3_trimm/ERR1698${i}_P1.fastq.gz -2 /sas/supercloud-kong/wangtianpeng/data/ath_flowering/3_trimm/ERR1698${i}_P2.fastq.gz -p 25 -o ./ERR1698${i}_quant & done

五、tximport的输入文件：涉及R的一些基本操作


dir<- getwd()
## "E:/Bio_informatics_software/转录组/NC文章重复_A.th"

list.files()
sample <- paste0("ERR1698", c(194, seq(202,209),seq(195,201)), "_quant")
### "ERR1698194_quant" "ERR1698202_quant" "ERR1698203_quant" "ERR1698204_quant"

files <- file.path(dir,sample,"quant.sf")
names(files) <- paste0("sample",1:16)
###                                                                               sample1 
"E:/Bio_informatics_software/转录组/NC文章重复_A.th/ERR1698194_quant/quant.sf" 
        sample2 
"E:/Bio_informatics_software/转录组/NC文章重复_A.th/ERR1698202_quant/quant.sf" 
        sample3 
"E:/Bio_informatics_software/转录组/NC文章重复_A.th/ERR1698203_quant/quant.sf"


file.exists(files)
all(file.exists(files))


library(AnnotationHub)
ah <- AnnotationHub()
### AnnotationHub with 40126 records
# snapshotDate(): 2017-04-25 
# $dataprovider: BroadInstitute, Ensembl
# $species:
# $rdataclass: GRanges, BigWigFile ...


ath <- query(ah, 'thaliana')
# snapshotDate(): 2017-04-25 
# $dataprovider:
# $species: 
# $rdataclass:
# additional mcols():

ath_tx <- ath[[ 'AH52247' ]]

k <- keys(ath_tx, keytype= "GENEID")
df <- select(ath_tx, keys= k, keytype= "GENEID", columns= "TXNAME")
tx2name <- df[,2:1]


library("tximport")
library("readr")
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)

# reading in files with read_tsv
# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 
# summarizing abundance
# summarizing counts
# summarizing length


names(txi)
### txi 为一个较大的list ，查看有哪些变量。"abundance"           "counts"              "length"  "countsFromAbundance
head(txi$counts)

六、DESeq2差异表达分析

dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition)

library(DESeq2)

### DESeqDataSetFromTximport(txi, colData, design)

condition<- factor(rep(c("DAY0", "DAY1", "DAY2", "DAY3"),each=4))
sampleTable <- data.frame(row.names = colnames(txi$counts),condition)
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~condition)

[sampleTable ]


image 

dim(dds) ### 显示32753个 16个sample
dds <- dds[rowSums(counts(dds))>1, ]
dim(dds)  ### 显示 24937个 16个sample

dds_out <- DESeq(dds)
res <- results(dds_out)
summary(res)



res0_1 <- results(dds_out, contrast = c("condition","DAY1","DAY0"))
res0_2 <- results(dds_out,contrast =c("condition","DAY2","DAY0"))
res0_3<- results(dds_out,contrast =c("condition","DAY3","DAY0"))

res0_1UP <- subset(res0_1,padj< 0.01 & log2FoldChange >0 )
res0_2UP <- subset(res0_2,padj< 0.01 & log2FoldChange >0 )
res0_3UP <- subset(res0_3,padj< 0.01 & log2FoldChange >0 )
res0_1DOWN <- subset(res0_1,padj< 0.01 & log2FoldChange <0 )
res0_2DOWN <- subset(res0_2,padj< 0.01 & log2FoldChange <0 )
res0_3DOWN <- subset(res0_3,padj< 0.01 & log2FoldChange <0 )

UP_GENE <- unique(rownames(res0_1UP),rownames(res0_2UP),rownames(res0_3UP))
DOWN_GENE <- unique(rownames(res0_1DOWN),rownames(res0_2DOWN),rownames(res0_3DOWN))
length(UP_GENE) ### 97个
length(DOWN_GENE)  ## 95个

七、 特殊基因的表达情况

library(AnnotationHub)
hub <- AnnotationHub()
query（hub, "thaliana"）
ath.tair.db <- ath[["AH53758"]]
keytypes(ath.tair.db)##查看有哪些数据类型，包含着各大主流数据库的数据。
###用select函数，就可以把任意公共数据库的数据进行一一对应。
### keys是原始的ID，columns是转换之后的ID，keytype是要指定的原始ID类型


plotGeneCounts <- function(genes) {
  require(ggplot2)
  require(tidyr)
  require(dplyr)
  GeneName <- select(ath.tair.db, keys=genes, columns=c("TAIR"),keytype = "SYMBOL")
  GeneName <- GeneName[order(GeneName$TAIR), ]
  Data <- dds[which(rownames(dds) %in% GeneName$TAIR)]
  Data1 <- as.data.frame(assay(Data))
  Data1$gene <- GeneName$SYMBOL[GeneName$TAIR == rownames(Data1)]
  Data2 <- Data1 %>%
    gather(sample,count, -gene) %>%
    mutate(condition = rep(c("Day0","Day1","Day2","Day3"), each= length(sample) / 4))
  p <- ggplot(Data2)
  p1 <- p + geom_boxplot(aes(x=gene,y=count, fill=condition))
  p2 <- p1 + theme(axis.title.x = element_blank()) + ylab("Counts")
  print(p2)
  return(GeneName)
}

-结果： 早期花的同源异型基因的表达量十分低。

上一篇 下一篇

猜你喜欢

热点阅读