RNAseq 教程 (2.2)

# 目录

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

# 2.2 Alignment

# HISAT2 alignment

用 HISAT2 比对基因组和转录组。

首先，为对齐结果创建适当的输出目录

mkdir align

HISAT2 的输出是每个数据集的 SAM/BAM 文件。

参考 HISAT2 帮助手册获得更多说明：

• https://ccb.jhu.edu/software/hisat2/manual.shtml

HISAT2 基本用法:

#hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <sam>]

额外参数如下:

• '-p 8' tells HISAT2 to use eight CPUs for bowtie alignments.

• '--rna-strandness RF' specifies strandness of RNAseq library. We will specify RF since the TruSeq strand-specific library was used to make these libraries. See here for options.

• '--rg-id $ID' specifies a read group ID that is a unique identifier.

• '--rg SM:$SAMPLE_NAME' specifies a read group sample name. This together with rg-id will allow you to determine which reads came from which sample in the merged bam later on.

• '--rg LB:$LIBRARY_NAME' specifies a read group library name. This together with rg-id will allow you to determine which reads came from which library in the merged bam later on.

• '--rg PL:ILLUMINA' specifies a read group sequencing platform.

• '--rg PU:$PLATFORM_UNIT' specifies a read group sequencing platform unit. Typically this consists of FLOWCELL-BARCODE.LANE

• '--dta' Reports alignments tailored for transcript assemblers.

• '-x /path/to/hisat2/index' The HISAT2 index filename prefix (minus the trailing .X.ht2) built earlier including splice sites and exons.

• '-1 /path/to/read1.fastq.gz' The read 1 FASTQ file, optionally gzip(.gz) or bzip2(.bz2) compressed.

• '-2 /path/to/read2.fastq.gz' The read 2 FASTQ file, optionally gzip(.gz) or bzip2(.bz2) compressed.

• '-S /path/to/output.sam' The output SAM format text file of alignments.

hisat2 -p 8 --rg-id=UHR_Rep1 --rg SM:UHR --rg LB:UHR_Rep1_ERCC-Mix1 --rg PL:ILLUMINA --rg PU:CXX1234-ACTGAC.1 -x INDEX/index --dta --rna-strandness RF -1 UHR_Rep1_ERCC-Mix1_Build37-ErccTranscripts-chr22.read1.fastq.gz -2 UHR_Rep1_ERCC-Mix1_Build37-ErccTranscripts-chr22.read2.fastq.gz -S align/UHR_Rep1.sam

  hisat2 -p 8 --rg-id=UHR_Rep2 --rg SM:UHR --rg LB:UHR_Rep2_ERCC-Mix1 --rg PL:ILLUMINA --rg PU:CXX1234-TGACAC.1 -x INDEX/index --dta --rna-strandness RF -1 UHR_Rep2_ERCC-Mix1_Build37-ErccTranscripts-chr22.read1.fastq.gz -2 UHR_Rep2_ERCC-Mix1_Build37-ErccTranscripts-chr22.read2.fastq.gz -S align/UHR_Rep2.sam

  hisat2 -p 8 --rg-id=UHR_Rep3 --rg SM:UHR --rg LB:UHR_Rep3_ERCC-Mix1 --rg PL:ILLUMINA --rg PU:CXX1234-CTGACA.1 -x INDEX/index --dta --rna-strandness RF -1 UHR_Rep3_ERCC-Mix1_Build37-ErccTranscripts-chr22.read1.fastq.gz -2 UHR_Rep3_ERCC-Mix1_Build37-ErccTranscripts-chr22.read2.fastq.gz -S align/UHR_Rep3.sam

  hisat2 -p 8 --rg-id=HBR_Rep1 --rg SM:HBR --rg LB:HBR_Rep1_ERCC-Mix2 --rg PL:ILLUMINA --rg PU:CXX1234-TGACAC.1 -x INDEX/index --dta --rna-strandness RF -1 HBR_Rep1_ERCC-Mix2_Build37-ErccTranscripts-chr22.read1.fastq.gz -2 HBR_Rep1_ERCC-Mix2_Build37-ErccTranscripts-chr22.read2.fastq.gz -S align/HBR_Rep1.sam

  hisat2 -p 8 --rg-id=HBR_Rep2 --rg SM:HBR --rg LB:HBR_Rep2_ERCC-Mix2 --rg PL:ILLUMINA --rg PU:CXX1234-GACACT.1 -x INDEX/index --dta --rna-strandness RF -1 HBR_Rep2_ERCC-Mix2_Build37-ErccTranscripts-chr22.read1.fastq.gz -2 HBR_Rep2_ERCC-Mix2_Build37-ErccTranscripts-chr22.read2.fastq.gz -S align/HBR_Rep2.sam

  hisat2 -p 8 --rg-id=HBR_Rep3 --rg SM:HBR --rg LB:HBR_Rep3_ERCC-Mix2 --rg PL:ILLUMINA --rg PU:CXX1234-ACACTG.1 -x INDEX/index --dta --rna-strandness RF -1 HBR_Rep3_ERCC-Mix2_Build37-ErccTranscripts-chr22.read1.fastq.gz -2 HBR_Rep3_ERCC-Mix2_Build37-ErccTranscripts-chr22.read2.fastq.gz -S align/HBR_Rep3.sam

注意：在上面的对齐中，我们将每个库视为一个独立的数据集。如果你有一个库的多个数据，你可以在一个 HISAT2 命令中将它们对齐在一起。要组合多个 reads 文件，您需要为 '-1' 输入参数提供所有 read1 文件作为逗号分隔的列表，然后为 '-2' 输入参数提供所有 read2 文件作为逗号分隔的列表，(其中两个列表的顺序相同): 还可以使用 samtool

# SAM 转换为 BAM

将 HISAT2 sam 文件转换为 bam 文件，并按对齐位置排序

# 合并 BAM 文件

将所有 UHR 数据和所有 HBR 数据合并成一个 BAM 文件。注意：这可以通过几种方式来完成，比如 “samtools merge”，“bamtools merge”，或者使用 picard-tools (见下文)。我们选择第三种方法是因为它在合并 bam 头信息方面做得最好。注意:sambamba 也保留头部信息。

cd align

java -Xmx2g -jar ../../picard.jar MergeSamFiles OUTPUT=UHR.bam INPUT=UHR_Rep1.bam INPUT=UHR_Rep2.bam INPUT=UHR_Rep3.bam

java -Xmx2g -jar ../../picard.jar MergeSamFiles OUTPUT=HBR.bam INPUT=HBR_Rep1.bam INPUT=HBR_Rep2.bam INPUT=HBR_Rep3.bam

计算对齐 (BAM) 文件，确保所有文件都成功创建 (总共应该有 8 个)

ls -l *.bam | wc -l

8

ls *.bam

HBR.bam       HBR_Rep2.bam  UHR.bam       UHR_Rep2.bam

HBR_Rep1.bam  HBR_Rep3.bam  UHR_Rep1.bam  UHR_Rep3.bam

# 练习 6

任务：对额外的数据集执行一些比对。用你在上面学到的技巧来对齐阅读。尝试使用 HISAT2。还要练习将 SAM 转换为 BAM 文件，并合并 BAM 文件。

在练习 3 中创建的名为 “practice” 的单独工作目录中进行分析。

.sam 和.bam 文件之间有什么区别？

sam 文件是一个纯文本序列比对映射文件。bam 文件是相同信息的二进制压缩版本。

如果您像上面所做的那样对结果 BAM 文件进行排序，那么结果是否按 read 名称排序？还是 position?

按照 position 排序

可以查看 BAM 文件的哪些列以确定排序的样式？

第一、第三和第四列包含 reads 名称、染色体和位置。

samtools view HCC1395_normal.bam | head | cut -f 1,3,4

可以使用什么命令仅查看 BAM 头？

samtools view -H HCC1395_normal.bam