以前寫過不少零散的 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 重復率會比較高慰照。ATCG 堿基占比容易出現(xiàn)前面堿基不穩(wěn)定。多查看幾個樣本表現(xiàn)琉朽，如果波動模式差不多說明是隨機引物不隨機導致的毒租，可以選擇移除或不移除。
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 文庫應有 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 對應哪個屬性萧福，比如 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è)置小點秘蛇∑渥觯或者比對率偏低時顶考，可以將 -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)含子有著復雜的關(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)有的基因重復了薄声，這是因為軟件認為有重復的基因当船。但提供的 GTF 這個基因并沒有重復，但是它存在不重疊的轉(zhuǎn)錄本默辨，軟件就認定為重復的基因德频，這算是軟件的一個 bug 了，希望哪天能修復了缩幸。
$ 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 導入給 DESeq2. 如果用 featureCounts 定量表谊，將 counts 整理成矩陣直接導入 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ù)導入到 DESeq2 分析媳板。
導入必要的 R 包和文件路徑桑腮。biomaRt 用于基因名轉(zhuǎn)換，因為差異分析我習慣用 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 文件路徑命名，保證讀取后矩陣的列跟文件是對應的提陶。
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 不應該進行過濾一些基因烫沙，除非在所有樣本 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)
導入到 DESeq2 并過濾低表達基因隙笆。過濾沒有統(tǒng)一標準锌蓄，我習慣要求至少在 n 樣本 counts 不小于 x. 如果數(shù)據(jù)有不同批次，一般在 design 公式將批次列（因子）放前面撑柔，注意 DESeq2 不會進行批次效應移除瘸爽，但分析時會區(qū)分哪些差異由批次效應引起，哪些由實驗條件引起乏冀。
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 測序至少有三個生物學重復辆沦，如果重復少于三個無法計算 cook's 距離。DESeq2 不支持技術(shù)重復识虚，如果做了技術(shù)重復肢扯，用 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 差異是否是實驗條件導致多糠，從而保留應有的差異避免過度壓縮累舷，得到的矩陣適合進行下游分析。
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é)論是相同的搭伤≈辉酰“壓縮”后的差異表達數(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 有詳細說明。
簡單來說肠骆，如果認為因子間有相互作用算途，不是獨立的。應在 design 參數(shù)添加相互作用項蚀腿。
假設(shè)實驗設(shè)計有 2 因子 genotype 值分別為 I/II/III 和 condition 值分別為 A/B.
如果 design 為 ~ genotype + condition 那么結(jié)果將是總體的 condition 導致的差異基因嘴瓤，排除了 genotype 因子的效應；如果是 ~ genotype + condition + genotype:condition 或者寫為 ~ genotype * condition 莉钙，那么 condition 因子效應將拆分為主要的效應——即 condition 因子在 genotype 為對照組（I）的情況下的效應廓脆，及相互作用項genotypeII.conditionB 和 genotypeIII.conditionB 表示對比于 genotypeI 在 genotypeII/genotypeIII 條件下的 condition 效應的差異。
下面舉個例子磁玉∧基因 1 在三種 genotype 下不同 condition 導致的差異倍數(shù)幾乎相同，那么相互作用項 genotype:condition 將接近于 0. 基因 2 在三種 genotype 下有不同的 condition 效應蜀涨，那么 genotype 為 II/III 的 condition 效應將是主要的 condition 效應加上對應 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
