RNA seq分析流程

2021-10-24 23:30:38 字數 3043 閱讀 8582

高通量測序知識

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 基本分析流程

easoncheng 高通量測序技術,就是二代測序,已經成為現代生物學研究的乙個較為常規的實驗手段。這一技術的發展極大地推動了基因組學,表觀基因組學以及翻譯組學的研究。rna seq 通過測定穩定狀態下的rna樣品的序列來對rna樣品進行研究,從而避免了許多之前研究手段的不足,比如象基因晶元或者 p...

RNA seq的典型流程(protocol)

一 rna的分離 從新鮮的或者是冷凍的細胞或組織樣本中分離rna,一般情況下,樣品會被dna汙染。因此,在製備文庫之前,會使用dna酶 降解dna 降解rna樣品中的dna汙染。二 rna的質量檢查 在後續分析的時候也有質控這一步,不過是從測序質量這個層面的質量 在製備文庫之前,要在rna降解,純度...

RNA Seq分析筆記(2)

檔案分割 使用下面命令將srr分割開來。fastq dump gzip split 3 srr3589956.sra done生成之後,結果如下 biochem qcgate step3 hisat2 ls gz l rw rw r 1 biochem biochem 1223462176 nov ...