RNAseq教程（3.4）

目录

1.Module 1 - Introduction to RNA sequencing

1. Installation
2. Reference Genomes
3. Annotations
4. Indexing
5. RNA-seq Data
6. Pre-Alignment QC

2.Module 2 - RNA-seq Alignment and Visualization

1. Adapter Trim
2. Alignment
3. IGV
4. Alignment Visualization
5. Alignment QC

3.Module 3 - Expression and Differential Expression

1. Expression
2. Differential Expression
3. DE Visualization
4. Kallisto for Reference-Free Abundance Estimation

4.Module 4 - Isoform Discovery and Alternative Expression

1. Reference Guided Transcript Assembly
2. de novo Transcript Assembly
3. Transcript Assembly Merge
4. Differential Splicing
5. Splicing Visualization

5.Module 5 - De novo transcript reconstruction

1. De novo RNA-Seq Assembly and Analysis Using Trinity

6.Module 6 - Functional Annotation of Transcripts

1. Functional Annotation of Assembled Transcripts Using Trinotate

3.4 Alignment Free Expression Estimation (Kallisto)

获取转录本的fasta序列

请注意，我们已经在RNA-seq教程的前面有了参考基因组序列的fasta序列。然而，Kallisto直接作用于目标cDNA/转录本序列。记住，我们有22号染色体上基因的转录模型。这些转录模型是以GTF格式从Ensembl下载的。GTF包含组成每个转录本的外显子的坐标描述，但不包含转录本序列本身。所以目前我们没有Kallisto需要的转录序列。我们有很多地方可以得到这些转录序列。

为了将Kallisto结果与StringTie的表达式结果进行比较，我们将创建一个定制的Fasta文件，该文件对应于用于StringTie分析的文本。我们如何以Fasta格式获得这些转录序列?

我们可以下载完整的人类fasta转录本数据库并只提取出22号染色体上的基因。或者我们可以使用来自tophat的名为gtf_to_fasta的工具从我们的GTF文件中生成fasta序列。这种方法很方便，因为它还包括控制中的ERCC峰值序列，允许我们为这些特征生成Kallisto丰度估计值。

gtf_to_fasta chr22_with_ERCC92.gtf chr22_with_ERCC92.fa chr22_ERCC92_transcripts.fa

使用less查看文件chr22_ERCC92_transcripts.fa。请注意，该文件有混乱的文字记录名称。使用下面的hairball perl一行代码来整理每个fasta序列的标题行

cat chr22_ERCC92_transcripts.fa | perl -ne 'if ($_ =~/^\>\d+\s+\w+\s+(ERCC\S+)[\+\-]/){print ">$1\n"}elsif($_ =~ /\d+\s+(ENST\d+)/){print ">$1\n"}else{print $_}' > chr22_ERCC92_transcripts.clean.fa
wc -l chr22_ERCC92_transcripts*.fa
  126025 chr22_ERCC92_transcripts.clean.fa
  126025 chr22_ERCC92_transcripts.fa
  252050 总用量

使用less chr22_ERCC92_transcripts.clean.fa查看生成的“clean”文件。使用tail查看文件末尾chr22_ERCC92_transcripts.clean.fa。请注意，我们对22号染色体上的每个集合转录本都有一个fasta记录，我们对每个ERCC spike-in序列都有一个额外的fasta记录。

cat chr22_ERCC92_transcripts.clean.fa | grep ">" | perl -ne '$_ =~ s/\>//; print $_' | sort | uniq > transcript_id_list.txt

Build a Kallisto transcriptome index

请记住，Kallisto并不进行比对或使用参考基因组序列。相反，它执行伪比对以确定读段与目标(本例中为转录序列)的兼容性。然而，与Tophat或STAR等对齐算法类似，Kallisto需要一个索引来高效、快速地评估这种兼容性。

kallisto index --index=chr22_ERCC92_transcripts_kallisto_index ../chr22_ERCC92_transcripts.clean.fa

Generate abundance estimates for all samples using Kallisto

正如使用StringTie和HT-Seq所做的那样，使用Kallisto为每个演示样本生成转录本丰度。

kallisto quant --index=chr22_ERCC92_transcripts_kallisto_index --output-dir=UHR_Rep1_ERCC-Mix1 --threads=4 --plaintext UHR_Rep1_ERCC-Mix1_Build37-ErccTranscripts-chr22.read1.fastq.gz UHR_Rep1_ERCC-Mix1_Build37-ErccTranscripts-chr22.read2.fastq.gz
kallisto quant --index=chr22_ERCC92_transcripts_kallisto_index --output-dir=UHR_Rep2_ERCC-Mix1 --threads=4 --plaintext UHR_Rep2_ERCC-Mix1_Build37-ErccTranscripts-chr22.read1.fastq.gz UHR_Rep2_ERCC-Mix1_Build37-ErccTranscripts-chr22.read2.fastq.gz
kallisto quant --index=chr22_ERCC92_transcripts_kallisto_index --output-dir=UHR_Rep3_ERCC-Mix1 --threads=4 --plaintext UHR_Rep3_ERCC-Mix1_Build37-ErccTranscripts-chr22.read1.fastq.gz UHR_Rep3_ERCC-Mix1_Build37-ErccTranscripts-chr22.read2.fastq.gz

kallisto quant --index=chr22_ERCC92_transcripts_kallisto_index --output-dir=HBR_Rep1_ERCC-Mix2 --threads=4 --plaintext HBR_Rep1_ERCC-Mix2_Build37-ErccTranscripts-chr22.read1.fastq.gz HBR_Rep1_ERCC-Mix2_Build37-ErccTranscripts-chr22.read2.fastq.gz
kallisto quant --index=chr22_ERCC92_transcripts_kallisto_index --output-dir=HBR_Rep2_ERCC-Mix2 --threads=4 --plaintext HBR_Rep2_ERCC-Mix2_Build37-ErccTranscripts-chr22.read1.fastq.gz HBR_Rep2_ERCC-Mix2_Build37-ErccTranscripts-chr22.read2.fastq.gz
kallisto quant --index=chr22_ERCC92_transcripts_kallisto_index --output-dir=HBR_Rep3_ERCC-Mix2 --threads=4 --plaintext HBR_Rep3_ERCC-Mix2_Build37-ErccTranscripts-chr22.read1.fastq.gz HBR_Rep3_ERCC-Mix2_Build37-ErccTranscripts-chr22.read2.fastq.gz

创建一个TSV文件，其中包含所有六个样本的TPM丰度估计值。

paste */abundance.tsv | cut -f 1,2,5,10,15,20,25,30 > transcript_tpms_all_samples.tsv
ls -1 */abundance.tsv | perl -ne 'chomp $_; if ($_ =~ /(\S+)\/abundance\.tsv/){print "\t$1"}' | perl -ne 'print "target_id\tlength$_\n"' > header.tsv
cat header.tsv transcript_tpms_all_samples.tsv | grep -v "tpm" > transcript_tpms_all_samples.tsv2
mv transcript_tpms_all_samples.tsv2 transcript_tpms_all_samples.tsv
rm -f header.tsv

head transcript_tpms_all_samples.tsv
tail transcript_tpms_all_samples.tsv