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