1介绍

本文描述了已知基因差异表达分析所需的一些“上游”步骤，使用超高速对齐工作流程，从FASTQ文件到计数数据，而无需中间BAM文件。

我们遵循鲑鱼教程,R在后来的阶段中被具体化。我们假设Linux操作系统带有文件夹〜/ a / salmon /教程， 与大马哈鱼安装在~ / bin /鲑鱼。

教程<- "~/a/salmon/ TUTORIAL " salmon <- "~/bin/salmon"

2下载

本教程要求FASTQ文件表示(配对端)序列样本。我们从短读存档(SRA)中检索它们wget命令行工具。

# !/bin/bash DIR = ~ / /教程SRA = ftp://ftp.sra.ebi.ac。Uk /vol1/fastq for I in ' seq 28 40 ';mkdir -p ${DIR}/data/DRR0161${i};cd $ {DIR} /数据/ DRR0161 $ {};wget -bqc $ {SRA} / DRR016 DRR0161 ${我}/ DRR0161 ${我}_1.fastq.gz;wget -bqc $ {SRA} / DRR016 DRR0161 ${我}/ DRR0161 ${我}_2.fastq.gz;做cd DIR美元

我们将序列与参考转录组对齐，在这种情况下，来自Ensembl拟南芥。我们下载序列并索引它。

# !/ bin / bash dir =〜/ a / tutorial ensembl = ftp：//ftp.ensemblinomes.org/pub/plants/release-28 curl \ $ {ensembl} /fasta/arabidopesis_thaliana/cdna/arabidopsis_thaliana.tair10.28.cdna.all.fa.gz \ -o $ {dir} / stathal.fa.gz salmon index -t $ {dir} / stathal.fa.gz -i $ {dir} / ahthal_index

如果我们感兴趣的话，很容易理解R

库(Biostrings) readDNAStringSet(文件。路径(教程中,“athal.fa.gz”))

##一个DNAStringSet实例的长度为41392 ## seq names ## [1] 462 ATGTCCCTTCTGTTTCAACA…GAAGAGATGTCGGCATTAA ATMG00010.1 cdna:…## [2] 324 atgtttcagttcgctaagtt…CTTTCTTTTCAACTCGTAA ATMG00030.1 cdna:…## [3] 948 atgacaaagcgtgagtataa…CAGACATTATGTCGACTAG ATMG00040.1 cdna:…## [4] 396 atgtctggcgtttaccggac…AAGCGGCGGGTTCGGGTAA ATMG00050.1 cdna:…## [5] 542 ccaatttgggccaattcc…AAAGTCAAGTCAAGAATAA ATMG00060.1 cdna:…## ... ... ... ## [41388] 952 GCCATTCCACATTTTTTTAT...ATGAAAAGAAGTAGGAATT AT1G80970.2 cdna:... ## [41389] 1418 GAACATAGTTAGACCGTGTA...CTTTTGTTTGTTGAATCAC AT1G80980.1 cdna:... ## [41390] 690 ATGGATAGATTGAGTGACTT...TCGTGACGACTGCTGCTGA AT1G80990.1 cdna:... ## [41391] 276 ATGTTATGTGCGTGTTGCGT...TGCTTTTTCGGATGATTGA AT1G81001.1 cdna:... ## [41392] 739 AAAGAAACCTCTTCCTCTAA...ATTAGAATAACATTATTTT AT1G81020.1 cdna:...

我们还需要知道转录本是如何映射到基因上的，因为我们要做基因水平的分析。我们下载匹配文本数据库的注释文件(gff3)。

# !/bin/bash DIR=~/a/tutorial ENSEMBL=ftp://ftp.ensemblgenomes.org/pub/plants/release-28 curl \ ${ENSEMBL}/gtf/arabidopsis_thaliana/ arabidopsis_thaliana . tair10.28 .gff3.gz\ -o ${DIR}/ arabidopsis_thaliana . tair10.28 .gff3.gz

3.量化

教程Bash脚本调用大马哈鱼在每个样本上，使用配对结束的FASTQ文件。

# !/home/ubuntu/bin/ SALMON=/home/ubuntu/bin/ SALMON for fn in $DIR/data/DRR0161{25..40};-i athal_index -l A \ -1 ${fn}/${sample}_1.fastq.gz \ -2 ${fn}/${sample}_2.fastq.gz \ -p 8 -o quants/${sample}_quant done

一个R量化一个样本的版本是

量化=函数（Sample_file，输出，鲑鱼，doit = true）{sample = baseName（sample_file）samples = file.path（sample_file，paste0（样本，“_”，1：2，“.fastq.gz”））输出= file.path（输出，paste0（样本，“_quant”））args = c（“Quant”，“-i”，索引，“-l a”，“-1”，样本[1]，“-2“，样本[2]，”-p“，并行:: detectcores（），”-o“，输出）if（doit）system2（salmon，args）else {txt = paste（salmon，paste（args，崩溃=“”））消息（txt）不可见（0）}}

在这里，我们量化所有的样本(集dooit = true来实际执行计数步骤。

指数=文件。path(TUTORIAL， "athal_index") output = file。data_dir = file. path(TUTORIAL， "quants")path(TUTORIAL， "data") sample_files = normalizePath(dir(data_dir, full = TRUE)) args = list(salmon = salmon, output = output, doit = FALSE) Map(quantify, sample_files, MoreArgs = args)

4分析R

我们需要提供解读所对应的转录本和我们将要进行分析的基因之间的映射。我们导入GFF注释文件，并提取转录本及其双亲(基因)。

Library (rtracklayer) Library (tidyverse) file = file。path(TUTORIAL， "Arabidopsis_thaliana.TAIR10.29.gff3") gff = import(file) tx2gene = tibble(txid = gff$transcript_id, geneid = as.character(gff$Parent))

库(summarizedexperexperiment) files = dir(文件。path(TUTORIAL， "quants")， pattern = "。sf"， recursive = TRUE, full = TRUE) names(files) = sub(“_quant”，“”，basename(dirname(files)))) counts (files, type =“salmon”，tx2gene = tx2gene)计数

##用READ_TSV读取文件

## 12 3 4 5 6 7 8 9 10 11 12 13 14 15 16总结的数量##总结的数量##总结的长度

名称（计数）

## [1]“丰富”“计数”“长度”## [4]“Countsfromabundance”

library(summarizedexperexperiment) se = summarizedexperexperiment (counts[-4])