以前寫過不少零散的 RNA-Seq 分析文章，現(xiàn)在整理為流程，同時修改一些錯誤。
流程包含質(zhì)控、比對、定量、差異分析。
前處理
拿到原始 fastq 數(shù)據(jù)先進行前處理。前處理包含質(zhì)控、比對和定量。質(zhì)控采用 fastqc/fastp; 比對用 hisat2 或者不比對，用 salmon 直接定量；比對后用 featureCounts 進行 reads 定量，用 TPMCalculator 進行 TPM 定量。
fastp 質(zhì)控
推薦 fastp 質(zhì)控原始 fastq 數(shù)據(jù)，因為集成程度高包含了各項過濾同時速度很快。
fastp -i ${RawDir}/${Sample}_R1.fq.gz -o ${CleanDir}/${Sample}_R1.fq.gz \
-I ${RawDir}/${Sample}_R2.fq.gz -O ${CleanDir}/${Sample}_R2.fq.gz \
--compression 6 --report_title ${Sample} --json ${CleanDir}/${Sample}_fastp.json \
--html ${CleanDir}/${Sample}_fastp.html --detect_adapter_for_pe
一般來說 RNA-Seq 重復(fù)率會比較高。ATCG 堿基占比容易出現(xiàn)前面堿基不穩(wěn)定。多查看幾個樣本表現(xiàn)，如果波動模式差不多說明是隨機引物不隨機導(dǎo)致的，可以選擇移除或不移除。
hisat2 比對
比對推薦 hisat2 不推薦 bowtie2. 因為 bowtie2 對剪切不友好，即使提高 --maxins 參數(shù)允許更長的片段，比對率依舊比較低。
hisat2 --new-summary -k 4 --threads 4 --summary-file ${BamDir}/${Sample}_HISAT2.txt \
-x ${GRCh38} -1 ${CleanDir}/${Sample}_R1.fq.gz -2 ${CleanDir}/${Sample}_R2.fq.gz \
-S ${BamDir}/${Sample}.sam
# samtools 過濾
samtools view --threads 4 -F 256 -b ${BamDir}/${Sample}.sam | samtools sort \
--threads 4 -o ${BamDir}/${Sample}.bam -O BAM
samtools index ${BamDir}/${Sample}.bam
ENCODE 2016 指南認為 Long/Small RNA-Seq 文庫應(yīng)有 30M 以上比對 Reads/Mates.
hisat2 默認 -k 參數(shù)為 5(linear index) 或 10(graph index)，此時 MAPQ 數(shù)值將為 0 或 1 或 60. 60 為唯一比對（primary），0 為不比對，1 為多比對。
$ samtools view Test.bam | awk '{print $5}' | sort | uniq
0
1
60
在 samtools view 命令用 -F 256 移除次要比對，保留唯一比對。
$ samtools flag 256
0x100 256 SECONDARY
featureCounts 定量
featureCounts 定量基礎(chǔ)單元叫 feature (如外顯子 exon)，以 GTF 注釋文件為例，其第三列為 feature 類型。多個 feature 可聯(lián)合形成 meta-feature (如基因 gene) 并計算總聚的 reads 數(shù)目。由 -g 參數(shù)設(shè)置 meta-feature 對應(yīng)哪個屬性，比如 GTF 格式默認用 "gene_id" 標記。
featureCounts 默認不計算重疊于多個 (meta-)feature 的 reads，除非重疊的 feature 屬于同一個 meta-feature 那么只計算一次。重疊于多個 meta-feature 對 RNA-Seq 實驗建議不計算，因為無法得知該 reads 究竟來自哪個基因；對于 CHIP-Seq 建議計算。因此用 featureCounts 定量轉(zhuǎn)錄本是不合適的，因為同一基因的轉(zhuǎn)錄本有大量重合區(qū)域。
featureCounts 用 BAM 文件 "NH" 標簽獲取多比對信息。默認只計算唯一比對（primary），不計算多比對。用 -M 參數(shù)控制計算多比對；用 --primary 參數(shù)控制只計算唯一比對。
featureCounts 默認 reads 有一個堿基與 feature 重合就計算一次，可通過 --minOverlap, --fracOverlap 等參數(shù)設(shè)定閾值。同時 reads 上的任意 gap(insertions, deletions, exon-exon junctions or structural variants) 也算數(shù)。
featureCounts 輸入的 BAM 文件不要求排序。下面 -p 和 --countReadPairs 參數(shù)指定是雙端測序且統(tǒng)計 fragment 而不是 reads.
featureCounts -p --countReadPairs --primary -F GTF -t exon -g gene_id \
-a ${Annot} -o ${CountDir}/${Sample}.tsv ${BamDir}/${Sample}.bam
salmon 定量
salmon 是轉(zhuǎn)錄本定量軟件，使用轉(zhuǎn)錄本定量主要優(yōu)點一是更準確，比如樣本同一基因使用不同轉(zhuǎn)錄本此時基因長度并不相等，直接基因定量不準確；二是方便進行可變剪切分析。
salmon 可以從 fastq 直接定量也可以從比對好的 bam 定量，本文使用 fastq 定量。
第一次使用 salmon 前先建立索引，建索引前先準備 decoys 文件。decoys 指不屬于人基因組但是測序又往往會檢測到的序列，有些 decoys 序列跟一些轉(zhuǎn)錄本很相似，明確 decoys 序列有助于準確定量。
用 GitHub - COMBINE-lab/SalmonTools: Useful tools for working with Salmon output 倉庫的腳本可以制作 decoys 文件。
bash generateDecoyTranscriptome.sh -j 8 -a gencode.v35.annotation.gtf \
-g GRCh38.primary_assembly.genome.fa -t gencode.v35.transcripts.fa -o hsa_decoy
用 decoys 目錄里的 fa 文件建立索引。
salmon index -t hsa_decoy/gentrome.fa -i hsa_transcripts_index \
--decoys hsa_decoy/decoys.txt -k 31 --gencode
設(shè)置 -k 31 適合 75 bp 及更長 reads，如果 reads 更短需要設(shè)置小點?；蛘弑葘β势蜁r，可以將 -k 設(shè)小些，提高靈敏度。
定量前注意兩點，一是 salmon 要求 reads 無序；二是注意文庫類型。
如上圖所示，根據(jù) reads 方向?qū)?RNA-Seq 文庫分為不同類型。salmon LIBTYPE 參數(shù)包含三個部分，一是 reads 的相對方向，分別用 I/O/M 表示；二是文庫有無鏈特異性，分別用 S/U 表示；三是，假如文庫有鏈特異性，指明鏈方向，分別用 F/R 表示。如果第二部分是 U 就不需要設(shè)置第三部分。下面是這些字母代表的全稱。
I = inward
O = outward
M = matching
S = stranded
U = unstranded
F = read 1 (or single-end read) comes from the forward strand
R = read 1 (or single-end read) comes from the reverse strand
A = automatically determine
salmon 用 F/R 表示鏈方向，其他軟件可能用 F2R1/F1R2 形式。
下圖是與 TopHat 比較
注意命令 -l 參數(shù)放 -1/-2 之前。
salmon quant -l IU -i ${IndexPath} -1 ${CleanDir}/${Sample}_R1.fq.gz \
-2 ${CleanDir}/${Sample}_R2.fq.gz -p 4 -g ${GRCh38GTF} \
-o ${SalmonDir}/${Sample}
轉(zhuǎn)錄本定量保存在 quant.sf 文件。有 -g 參數(shù)時會輸出基因水平定量。
Name Length EffectiveLength TPM NumReads
ENST00000456328.2 1657 1372.105 0.034568 1.000
ENST00000450305.2 632 348.480 0.000000 0.000
ENST00000488147.1 1351 1066.105 3.604749 81.025
注意 NumReads 不一定是整數(shù)，是考慮唯一比對和多比對后對轉(zhuǎn)錄本相對豐度的估計。
TPMCalculator 定量
TPMCalculator 計算基因、轉(zhuǎn)錄本、外顯子和內(nèi)含子的 TPM 定量。其 TPM 定量公式如下， 是比對到
feature 片段數(shù)目；
是 feature 長度。
基因的轉(zhuǎn)錄本及外顯子內(nèi)含子有著復(fù)雜的關(guān)系，如下圖所示軟件進行了兩種轉(zhuǎn)換。第一種轉(zhuǎn)換，根據(jù)外顯子重合，創(chuàng)建包含所有外顯子的基因模型；第二種轉(zhuǎn)換創(chuàng)建“純粹”的內(nèi)含子區(qū)域，不被任何外顯子覆蓋。這些不重合的內(nèi)含子和 feature 用來計算 TPM.
參數(shù) -c 推薦設(shè)置為 reads 長度；參數(shù) -p 設(shè)置只統(tǒng)計“正確”比對的雙端測序數(shù)據(jù)。
因為比對 hisat2 設(shè)置了 -k 參數(shù)，因此 MAPQ 值只有 0/1/60 三種，用 -q 參數(shù)過濾比對質(zhì)量 60 因此只保留 primary alignment. 如果前面經(jīng)過 samtools 過濾了，這里不再需要過濾。
TPMCalculator -g ${GRCh38GTF} -b ${BamDir}/${Sample}.bam \
-c 150 -k gene_id -t transcript_id -q 60 -p -a
軟件將輸出結(jié)果至當前目錄，包含 3 文件。
- *_genes.out | Include calculated values per Gene
- *_genes.uni | Include calculated values per non-overlapped features per Gene
- *_genes.ent | Include calculated values for all features per Gene
仔細檢查基因輸出結(jié)果可能會發(fā)現(xiàn)有的基因重復(fù)了，這是因為軟件認為有重復(fù)的基因。但提供的 GTF 這個基因并沒有重復(fù)，但是它存在不重疊的轉(zhuǎn)錄本，軟件就認定為重復(fù)的基因，這算是軟件的一個 bug 了，希望哪天能修復(fù)了。
$ grep "ENSG00000235538" kLEC96h-2_genes.out | awk '{NF=7;print}'
ENSG00000235538.3#1 chr6 163671576 163798848 127273 4 0.00366198
ENSG00000235538.3#2 chr6 163927121 164009979 82859 2 0.00281244
ENSG00000235538.3#3 chr6 164228330 164231609 3280 0 0
差異表達分析
差異分析用 DESeq2 包，定量采用 salmon 結(jié)果，需要 tximport 將 salmon 導(dǎo)入給 DESeq2. 如果用 featureCounts 定量，將 counts 整理成矩陣直接導(dǎo)入 DESeq2 分析。
tximport 將 salmon 定量的轉(zhuǎn)錄本表達計算成基因表達。這需要基因和轉(zhuǎn)錄本映射表（TxDb），TxDb 用 GenomicFeatures 包從 GTF 注釋文件創(chuàng)建，然后保存到本地，以后就讀取使用。
library(AnnotationDbi)
library(GenomicFeatures)
gtfPath <- file.path("/example_dir/gencode.v35.annotation.gtf")
metaName <- c("Source", "Version", "Species")
metaValue <- c("GENCODE", "v35", "Homo sapiens")
extraInfo <- data.frame(name=metaName, value=metaValue)
txd <- makeTxDbFromGFF(gtfPath, format = "gtf", dataSource = "gencode.v35.annotation.gtf", organism = "Homo sapiens", metadata = extraInfo)
# 保存 TxDb 對象到本地
txdbPath <- file.path("/example_dir/gencode.v35.annotation.TxDb.sqlite")
saveDb(txd, file = txdbPath)
用 makeTxDbFromGFF 讀取 GTF 注釋到 TxDb. 參數(shù) dataSource 記載數(shù)據(jù)來源，可填可不填；參數(shù) metadata 提供元信息，格式是兩列的數(shù)據(jù)框，要求列名分別為 "name", "value".
> keytypes(txd)
[1] "CDSID" "CDSNAME" "EXONID" "EXONNAME" "GENEID" "TXID" "TXNAME"
> txToGene <- AnnotationDbi::select(txd, keys=keys(txd, "TXNAME"), keytype="TXNAME", columns=c("TXNAME", "GENEID"))
'select()' returned 1:1 mapping between keys and columns
> head(txToGene)
TXNAME GENEID
1 ENST00000456328.2 ENSG00000223972.5
2 ENST00000450305.2 ENSG00000223972.5
3 ENST00000473358.1 ENSG00000243485.5
4 ENST00000469289.1 ENSG00000243485.5
5 ENST00000607096.1 ENSG00000284332.1
6 ENST00000606857.1 ENSG00000268020.3
txToGenePath <- file.path("/example_dir/gencode.v35.annotation.TxToGene.csv")
write.csv(txToGene, file = txToGenePath, quote = FALSE, row.names = FALSE)
選擇轉(zhuǎn)錄本和基因兩列，得到他們映射關(guān)系，并保存到本地 csv 文件，下次直接讀取使用。注意第一列為轉(zhuǎn)錄本 ID 第二列為基因 ID 固定順序，列名不作要求。
有 TxDb 和 salmon 定量文件和樣品分組信息，就可以將數(shù)據(jù)導(dǎo)入到 DESeq2 分析。
導(dǎo)入必要的 R 包和文件路徑。biomaRt 用于基因名轉(zhuǎn)換，因為差異分析我習(xí)慣用 Ensembl ID 下游分析又常用 HGNC SYMBOL 或 ENTREZ ID.
library(tidyverse)
libraty(biomaRt)
library(tximport)
library(DESeq2)
groupPath <- file.path("/example_dir/SampleGroup.csv")
txGenePath <- file.path("/example_dir/gencode.v35.annotation.TxToGene.csv")
salmonDir <- file.path("/example_dir/Salmon")
ensembl <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", dataset = "hsapiens_gene_ensembl")
tximport 讀取 salmon 定量，并保存一份 TPM 表達矩陣，以備下游分析需要。要注意給 salmon 文件路徑命名，保證讀取后矩陣的列跟文件是對應(yīng)的。
sampleGroup <- read.csv(groupPath, header = TRUE, quote = "", row.names = 1)
sampleList <- rownames(sampleGroup)
sampleFiles <- file.path(salmonDir, sampleList, "quant.sf")
# 注意給路徑命名
# 名字是表達矩陣樣品名
names(sampleFiles) <- sampleList
txToGene <- read.csv(txGenePath, header = TRUE, quote = "", stringsAsFactors = FALSE)
txi <- tximport(sampleFiles, type = "salmon", tx2gene = txToGene)
# 理論上 TPM 不應(yīng)該進行過濾一些基因，除非在所有樣本 TPM 都為 0
tpm1 <- txi$abundance
tpm2 <- tpm1[rowSums(tpm1) > 0, ]
# 去除 Ensembl ID 版本號
short_id <- function(long_id) {
id_version <- strsplit(long_id, split = ".", fixed = TRUE) %>% unlist()
ensembl_id <- id_version[1]
return(ensembl_id)
}
# 取得所有表達基因，創(chuàng)建基因 ID 映射表
ensemblGenes <- sapply(rownames(tpm2), short_id, USE.NAMES = FALSE)
geneMap <- getBM(mart = ensembl, attributes = c("ensembl_gene_id", "entrezgene_id", "hgnc_symbol"),
values = ensemblGenes, filters = "ensembl_gene_id") %>%
as_tibble() %>%
arrange(entrezgene_id, desc(hgnc_symbol)) %>%
distinct(ensembl_gene_id, .keep_all = TRUE)
# 保存 TPM 矩陣到文件
tpm3 <- as_tibble(tpm2, rownames="gene_id") %>%
tidyr::separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
dplyr::select(-ensembl_gene_version) %>% dplyr::left_join(geneMap, by="ensembl_gene_id") %>%
dplyr::select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
tpmPath <- file.path("/example_dir/TPM.csv")
write_csv(tpm3, tpmPath)
導(dǎo)入到 DESeq2 并過濾低表達基因。過濾沒有統(tǒng)一標準，我習(xí)慣要求至少在 n 樣本 counts 不小于 x. 如果數(shù)據(jù)有不同批次，一般在 design 公式將批次列（因子）放前面，注意 DESeq2 不會進行批次效應(yīng)移除，但分析時會區(qū)分哪些差異由批次效應(yīng)引起，哪些由實驗條件引起。
dds1 <- DESeqDataSetFromTximport(txi, colData = sampleGroup, design = ~ Group)
# 過濾低表達基因
# 一共 6 樣品，要求至少 2 樣品 read counts 大于等于 5
keepGene <- rowSums(counts(dds1) >= 5) >= 2
dds2 <- dds1[keepGene, ]
dds3 <- DESeq(dds2)
用 cook's 距離檢測組內(nèi)樣本一致性。一般 RNA-Seq 測序至少有三個生物學(xué)重復(fù)，如果重復(fù)少于三個無法計算 cook's 距離。DESeq2 不支持技術(shù)重復(fù)，如果做了技術(shù)重復(fù)，用 collapseReplicates 函數(shù)合并。
ddsAssay <- assays(dds3)
cooks <- ddsAssay$cooks
boxplot(log10(cooks))
取得差異分析結(jié)果前保存表達矩陣，以備數(shù)據(jù)質(zhì)控和下游分析。用 counts 函數(shù)取得原始的 counts 或標準化的 counts. 后面差異分析結(jié)果的 baseMean 列就是標準化的 counts 計算得到。函數(shù) rlog 或 vst 將 counts 根據(jù)文庫大小或其他標準化因子進行了標準化和轉(zhuǎn)換到對數(shù)（log2）尺度，并控制方差獨立于表達均值。參數(shù) blind 控制轉(zhuǎn)換過程是否考慮實驗設(shè)計——即前面的 design 參數(shù)。設(shè)置 blind = FALSE 會考慮 counts 差異是否是實驗條件導(dǎo)致，從而保留應(yīng)有的差異避免過度壓縮，得到的矩陣適合進行下游分析。
readCounts <- DESeq2::counts(dds3, normalized = FALSE) %>%
as_tibble(rownames = "gene_id") %>%
separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
select(-ensembl_gene_version) %>%
left_join(geneMap, by = "ensembl_gene_id") %>%
select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
normalizedCounts <- DESeq2::counts(dds3, normalized = TRUE) %>%
as_tibble(rownames = "gene_id") %>%
separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
select(-ensembl_gene_version) %>%
left_join(geneMap, by = "ensembl_gene_id") %>%
select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
# 適合質(zhì)控
rLog1 <- rlog(dds3, blind = TRUE) %>% assay() %>%
as_tibble(rownames = "gene_id") %>%
separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
select(-ensembl_gene_version) %>%
left_join(geneMap, by = "ensembl_gene_id") %>%
select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
# 適合下游分析
rLog2 <- rlog(dds3, blind = FALSE) %>% assay() %>%
as_tibble(rownames = "gene_id") %>%
separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
select(-ensembl_gene_version) %>%
left_join(geneMap, by = "ensembl_gene_id") %>%
select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
取得差異分析結(jié)果并保存 MA 圖。其中 results 函數(shù)取得“原始”的差異表達數(shù)據(jù)，lfcShrink 將差異倍數(shù)“壓縮”，兩個數(shù)據(jù) p 值相同——即統(tǒng)計檢驗的結(jié)論是相同的?！皦嚎s”后的差異表達數(shù)據(jù)適合用于排序或可視化。要注意如果實驗設(shè)計 design 包含了交互作用項，不要進行 lfcShrink 因為容易過度“壓縮”。
參數(shù) contrast 接受多種指定方式，下面代碼是常見的兩種方式。第一種是三個元素的向量，向量第一個元素是因子名字，第二個元素是差異分析分子項的水平，即常指的實驗組（比較組），第三個元素是差異分析分母項的水平。舉例來說，c("Condition", "B", "A") 表示條件 B 樣本對比條件 A 樣本的差異結(jié)果。第二種從 resultsNames() 結(jié)果選取一個字符串，如 "Group_48H_vs_0H" 表示 Group 因子項的 48H 樣本比較 0H 樣本結(jié)果。
rawRes <- results(dds3, contrast = c("Group", "48H", "0H"))
lfcRes <- lfcShrink(dds3, coef = "Group_48H_vs_0H", type = "apeglm", res = rawRes)
summary(rawRes)
# MA Plot
plotMA(rawRes, ylim = c(-5, 5), main = "Raw")
plotMA(lfcRes, ylim = c(-5, 5), main = "Shrink")
比較 MA 圖和火山圖，可以看出 lfcShrink 對低表達基因“壓縮”影響更大。
將差異分析結(jié)果轉(zhuǎn)換為 tibble 并保存，下游分析根據(jù)需要進行過濾。一般來說取 "abs(log2FoldChange) >= 1, padj < 0.05" 作為差異表達基因。
rawDEGs <- as_tibble(rawRes, rownames = "gene_id") %>%
separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
select(-ensembl_gene_version) %>%
left_join(geneMap, by="ensembl_gene_id") %>%
select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
lfcDEGs <- as_tibble(lfcRes, rownames = "gene_id") %>%
separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
select(-ensembl_gene_version) %>%
left_join(geneMap, by="ensembl_gene_id") %>%
select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
差異基因結(jié)果基因 p 值可能是 NA. 包含下列情況：
- 該基因所有樣本 counts 數(shù)都為 0, 那么 baseMean 列為 0, 此時 log2FC, p, padj 為 NA
- 該基因包含離群的 counts 數(shù)，此時 p, padj 為 NA
- 該基因 baseMean 非常低，此時 padj 為 NA.
相互作用項
上面提到實驗設(shè)計包含相互作用項時不建議 lfcShrink 處理，相互作用項在官網(wǎng) Interactions 有詳細說明。
簡單來說，如果認為因子間有相互作用，不是獨立的。應(yīng)在 design 參數(shù)添加相互作用項。
假設(shè)實驗設(shè)計有 2 因子 genotype 值分別為 I/II/III 和 condition 值分別為 A/B.
如果 design 為 ~ genotype + condition 那么結(jié)果將是總體的 condition 導(dǎo)致的差異基因，排除了 genotype 因子的效應(yīng)；如果是 ~ genotype + condition + genotype:condition 或者寫為 ~ genotype * condition ，那么 condition 因子效應(yīng)將拆分為主要的效應(yīng)——即 condition 因子在 genotype 為對照組（I）的情況下的效應(yīng)，及相互作用項genotypeII.conditionB 和 genotypeIII.conditionB 表示對比于 genotypeI 在 genotypeII/genotypeIII 條件下的 condition 效應(yīng)的差異。
下面舉個例子?；?1 在三種 genotype 下不同 condition 導(dǎo)致的差異倍數(shù)幾乎相同，那么相互作用項 genotype:condition 將接近于 0. 基因 2 在三種 genotype 下有不同的 condition 效應(yīng)，那么 genotype 為 II/III 的 condition 效應(yīng)將是主要的 condition 效應(yīng)加上對應(yīng) genotype 的相互作用項。
參考資料
Vera Alvarez, Roberto, et al. "TPMCalculator: one-step software to quantify mRNA abundance of genomic features." Bioinformatics 35.11 (2019): 1960-1962.
QC Fail Sequencing ? MAPQ values are really useful but their implementation is a mess
Analyzing RNA-seq data with DESeq2
Rsubread/Subread Users Guide
Overview – Salmon: Fast, accurate and bias-aware transcript quantification from RNA-seq data
Bowtie 2: fast and sensitive read alignment
The decoy genome
https://www.biostars.org/p/456231/
Decoy Sequences V/S Target Sequences
Importing transcript abundance with tximport
