1.fastqc
2.star
##build_index
star --runthreadn 9 --runmode genomegenerate \
--genomedir /data/***xx/bio/task_le-mir/03-2mirseq/index \
--genomefastafiles /data/***xxbio/task_le-mir/03-2mirseq/chrom.37.fa \
--sjdbgtffile /data/***xx/bio/task_le-mir/03-2mirseq/mirna37note.gtf \
--sjdboverhang 149
###star_align
ls -d le*|while read le; do echo $le; star --runthreadn 5 --genomedir /data/***xx/bio/task_le-mir/03humansequence/index \
--readfilescommand zcat \
--readfilesin /data/***xx/bio/task_le-mir/01trim_out/$/*_1.fq.gz \
/data/***xx/bio/task_le-mir/01trim_out/$/*_2.fq.gz \
--outfilenameprefix /data/***xx/bio/task_le-mir/04align_out/$_ \
--outsamtype bam sortedbycoordinate \
--outbamsortingthreadn 5 \
--quantmode transcriptomesam genecounts>>$.log; done
###rsem-prepare
rsem-prepare-reference --gtf /data/***xx/bio/task_le-mir/03-2mirseq/mirna37note.gtf \
/data/***xx/bio/task_le-mir/03-2mirseq/chrom.37.fa \
/data/***xx/bio/task_le-mir/05-2rsem/rsem_prepare
##rsem calculate expression
ls -d le*|while read le; do echo $le; rsem-calculate-expression --paired-end --no-bam-output \
--alignments -p 10 \
-q $/*_aligned.totranscriptome.out.bam \
/data/***xx/bio/task_le-mir/05-2rsem/rsem_prepare \
/data/***xx/bio/task_le-mir/05-2rsem/rsem_out/$_ >>$.log; done
####將表達定量結果轉換為矩陣
rsem-generate-data-matrix le*_.isoforms.results >output.matrix
###deseq2.r
setwd("/data/***xx/bio/task_le-mir/06deseq2")
##讀取檔案
input_data <- read.table("deseq2_input-2.txt", header=true, row.names=1)
##取整
input_data <-round(input_data, digits = 0)
##1-i型 2-ii型 3-正常+良性 4-補充良性
1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 1, 1, 1, 1, 2, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 1, 1, 1, 1, 1, 2,
3, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 2, 3, 2, 2, 1, 1, 3, 3,
1, 4, 4, 4, 4, 4, 4, 4
##準備工作
input_data <- as.matrix(input_data)
condition <- factor(c(1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 1, 1, 1, 1, 2, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 2,
2, 3, 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1,
2, 3, 2, 2, 1, 1, 3, 3, 1, 4, 4, 4, 4, 4, 4, 4))
coldata <- data.frame(row.names=colnames(input_data), condition)
##載入包
library(deseq2)
##構建dds矩陣
dds <- deseqdatasetfrommatrix(countdata=input_data, coldata=coldata, design=~condition)
##差異分析
dds <- deseq(dds)
##提取結果
res <- results(dds)
##看結果
summary(res)
##按p值排序
res <- res[order(res$padj), ]
resdata <-merge(as.data.frame(res), as.data.frame(counts(dds,normalized=true)), by="row.names", sort=false)
names(resdata)[1] <- "isoform"
##輸出結果檔案
write.table(resdata, file="diffexpr-results.txt", sep="\t", quote=f, row.names=f)
##視覺化展示
plotma(res)
##提取差異結果
awk '' diffexpr-i,ii.txt > upgene.txt
awk '' diffexpr-i,ii.txt > downgene.txt
轉錄組分析 轉錄組分析 使用STAR進行比對
通過二代測序我們可以獲得150bp左右的reads，如果想要知道reads是從哪個轉錄本上測出來的，就需要將reads比對到參考基因組上。比對的演算法很複雜，但簡單理解就是看reads與基因組上哪個區域一致。wget c 解壓 star tar xvzf 2.7.3a.tar.gz 執行 star ...
轉錄組分析 高階轉錄組分析和R資料視覺化
常規轉錄組是我們最常接觸到的一種高通量測序資料型別，其實驗方法成熟，花費較低，是大部分cns必備的技術，以後應該就如做個pcr一樣常見。而且分析思路簡潔清晰，是入門生信，學習生信分析思路和資料視覺化的首選。資料分析是相通的，通過乙個簡單的課程理解其中的原理，就可以推而廣之，延伸到其它型別的資料分析，...
轉錄組分析的正確姿勢
轉錄組分析是目前應用最廣的高通量測序分析技術之一。常見設計是不同樣品之間比較，尋找差異基因 標誌基因 協同變化基因 差異剪接和新轉錄本，並進行結果視覺化 功能注釋和網路分析等。轉錄組的測序分析也相對成熟，從rna提取 構建文庫 上機測序再到結果解析既可以自己完成，又可以在專業公司進行。概括來看轉錄組...