1. 對(duì)gtf文件進(jìn)行處理：
cat GSE52313_transcripts.gtf | grep "gene_name" |sed 's/\(.*\)gene_id\(.*\); transcript_id\(.*\); gene_name\(.*\); oId\(.*\)/\2\4/g' |sed 's/"http://g' | awk '{print $1"="$2}' | sort -n | uniq -c > GSE52313_transcripts_sort.gtf
awk -F" *" '{print $3}' GSE52313_transcripts_sort.gtf | awk -F"=" '{print $1,$2}' > GSE52313_transcripts_renew.gtf
2. 提取lncRNA:
awk 'NR==FNR{a[$1]=$2}NR!=FNR{if ($1 in a){print a[$1],$2,$3,$4,$5,$6,$7,$8,$9}}' GSE52313_transcripts_renew.gtf GSE52313_lnRNA
3. 對(duì)lncRNA表達(dá)量文件進(jìn)行基因名替換
awk 'NR==FNR{a[$1]=$2}NR!=FNR{if ($1 in a){print a[$1],$2,$3,$4,$5,$6,$7,$8,$9}}' GSE52313_transcripts_renew.gtf GSE52313_lnRNA_gene_counts > GSE52313_lnRNA_gene_counts_renew
4.根據(jù)table S1的信息：sharm:?8R 3RL 12RL 7LR
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? MI:5O 7R 11L 2L
5. 對(duì)lncRNA進(jìn)行差異表達(dá)分析：
首先需要讀入一個(gè)數(shù)據(jù)框搂漠，列代表每個(gè)sample承绸，行代表每個(gè)gene
database_all <- read.table(file = "GSE52313_lnRNA_gene_counts_renew", sep = "\t", header = T, row.names=1)
這里主要對(duì)于兩兩比較的數(shù)據(jù)钠怯，因此我取了數(shù)據(jù)的前6列，分別是兩組樣品冕末，每組3個(gè)生物學(xué)重復(fù)
設(shè)定分組信息，也就是樣本分組的名稱(chēng)
type <- factor(c(rep("LC_1",4), rep("LC_2",4)))
我這里是樣品1是LC_1，樣品2是LC_2
由于DESeq包要求接下來(lái)的count data必須要整數(shù)型饮醇，因此我們需要對(duì)數(shù)據(jù)進(jìn)行取整,然后將數(shù)據(jù)database和分組信息type讀入到cds對(duì)象中
database <- round(as.matrix(database_all))
source("https://bioconductor.org/biocLite.R")
biocLite('BiocInstaller')
biocLite("DESeq2")
coldata <- data.frame(row.names = colnames(database), type)
dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~type)
DESeq2的使用方法：
輸入矩陣數(shù)據(jù)，行名為sample浪漠，列名為gene陕习；DESeq2不支持無(wú)生物學(xué)重復(fù)的數(shù)據(jù)，因此我選擇了2個(gè)樣本郑藏，3個(gè)生物學(xué)重復(fù)的數(shù)據(jù)衡查；并對(duì)count data取整（經(jīng)大神指點(diǎn)，這里需要說(shuō)明下必盖，我的測(cè)試數(shù)據(jù)readcount是RSEM定量的結(jié)果拌牲，并不是常見(jiàn)的htseq-count的結(jié)果，所以count值會(huì)有小數(shù)點(diǎn)歌粥，而DESeq2包不支持count數(shù)有小數(shù)點(diǎn)塌忽，所以這里需要round取整）。
設(shè)置分組信息以及構(gòu)建dds對(duì)象
dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)
使用DESeq函數(shù)進(jìn)行估計(jì)離散度失驶，然后進(jìn)行標(biāo)準(zhǔn)的差異表達(dá)分析土居，得到res對(duì)象結(jié)果
前過(guò)濾:
keep<-rowSums(counts(dds)) >= 10
dds<-dds[keep,]
注：并不是必須，不影響計(jì)算結(jié)果嬉探，這樣做的優(yōu)點(diǎn)是dds對(duì)象占用的內(nèi)存小點(diǎn)擦耀，后續(xù)的計(jì)算耗時(shí)小點(diǎn)。
dds <- DESeq(dds)
res <- results(dds)
最后設(shè)定閾值涩堤，篩選差異基因眷蜓，導(dǎo)出數(shù)據(jù)
table(res$padj <0.05)
res <- res[order(res$padj),]
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by="row.names",sort=FALSE)
write.csv(resdata,file = "LC_1_vs_LC_2.csv")
#從LC_1_vs_LC_2.csv中篩選lncRNA和mRNA，
