RNA seq分析流程

高通量測序知識

fastqc使用，相對應的r包fastqcr，

rqc

fastp

biostrings包計算gc含量，q20等

library(biostrings)
filepath <- system.file("extdata", "s_1_sequence.txt",
package="biostrings")
qdna2 <- readqualityscaleddnastringset(filepath)
qdna2
#得到每一行的gc含量
gc_content <- letterfrequency(dnastringset(qdna2), letters="cg",as.prob = t)
#得到整個fastq的gc含量
letterfrequency(dnastringset(qdna2), letters="cg",as.prob = t,collapse = t)
#q20，將質量分數轉為數值去計算,整個檔案
qa <- as(quality(qdna2),"integerlist")
library(dplyr)
sum_20 <- purrr::map_int(seq_along(qa),~length(which(qa[[.x]]>20)) %>% 
sum()
sum_all <- quality(qdna2) %>% width() %>% sum()
q20 <- sum_20/sum_all

在比對檔案之前需要先要去掉adapter，得到不含有adapter 的fastq 檔案。有乙個flexbark可以，

flexbar -t 輸出檔名 -qf i1.8 -r 輸入檔案 -a adapters.fa

我沒試過去接頭。

反正有了trimmed 資料後，就要與參考檔案進行比對。目前大家所認同的比對演算法

有很多，比如bowtie2，bwa，tophat，hisat2 和star 等等。試下

綜合指標 (sahraeian 等，2017) 比較高的star 進行比對。

比對大部分需要建立索引,加快速度。**

star --runthreadn 執行緒數 --runmode genomegenerate --genomedir 存放參考基因的目錄 --genomefastafiles 參考基因組的 fasta 檔案 --sjdbgtffile 參考基因組的gtf 檔案 --sjdboverhang

索引長度,這個是reads長度的最大值減1，預設是100

建立好索引就可以對資料進行比對了，比對如下：

star --genomedir 存放參考基因的目錄 --readfilesin 輸入檔案(read1,if need;read2) --runthreadn 執行緒數 --outsamtype bam --outfilenameprefix 輸

出檔案字首

具體可以參考

bed檔案從gtf的獲取

bedtools multicov -bams *.bam -bed bed檔案 >reads_by_bed.txt

雖然-bams可以使用萬用字元，但是我一般還是用shell寫指令碼迴圈，因為輸出檔案沒有列名，很容易混淆的。我認為輸出檔案最好能帶個列名。

user guide

比如雙端測序，

直接使用gene_name做識別符號也行,提取出gene的位置。

類似的還有htseq-count

而samtools 的bedcov計算的則不同，是bed檔案每個區域的深度和。即samtools depth的和。

而計算覆蓋到整個參考fasta可以用samtools flagstat或者更詳細的samtools idxstats。雖然featurecounts等也會給出，不過是輸出在標準輸出裡的，不在輸出檔案裡。

reads數即counts數，接著如果需要比較不同樣本同個基因上的表達豐度情況，則需要對count數進行標準化。

落在乙個基因區域內的read counts數目一般可以認為取決於length of the gene（基因長度）和sequencing depth（測序深度）

所以，這個fragments是不是featurecounts算出的counts數我不太確定。反正rpkm肯定是。

在r中，則

counts_to_fpkm <- function(counts,lengths)

而tpm則是

counts_to_tpm <- function(counts,lengths)

那麼，exon length呢，粗暴地理解為gene length？【這樣方便獲取】，正常來說應該是每個gene id的exon的長度總和。獲取所有外顯子的長度和呢，一般可以從gtf檔案入手，但是gtf檔案中的exon的區域很多都是重合的（不同轉錄本的可變剪下）【比較麻煩。】

當然，使用featurecounts時，可以

輸出檔案有一列length，即為exon的length,方便。

而在r中,可以使用genomicfeatures計算每個基因下所有外顯子的總長度

#製作txdb物件
library(genomicfeatures)
txdb <- maketxdbfromgff("hg38.gtf",format="gtf")
#通過exonsby獲取每個gene上的所有外顯子的起始位點和終止位點，然後用reduce去除掉重疊冗餘的部分，最後計算長度
exons_gene <- exonsby(txdb, by = "gene")
exons_gene_lens <- purrr::map_int(exons_gene,~sum(width(reduce(.x))))

差異分析使用的是counts檔案，相關性分析之類的可以使用tpm資料。

RNA seq分析流程

RNA seq 基本分析流程

RNA seq的典型流程（protocol）

RNA Seq分析筆記（2）

RNA seq分析流程

RNA seq 基本分析流程

RNA seq的典型流程（protocol）

RNA Seq分析筆記（2）

相關推薦