2.Expression Matrix

Pipeline

Data Structure

"genomes/hg38/"         # reference genomes (i.e. genome sequence and annotation)
"shared_scripts/"       # shared scripts by Lu Lab
~/github                # I sync my own scripts to github

~/proj_exRNA/
|-- RNA_index           # rerference transcriptomes (fasta and index) 
|-- sample*_name        # ...
|-- sample2_name        # different samples     
`-- sample1_name        
    |-- fastq           # raw data: fastq files
    |-- fastqc          # QC of fastq
    |-- trim            # trimmed fastq (e.g. 3' adaptor cutted)
    |-- mapped          # mapped data: SAM/BAM files
    `-- exp             # expression of each gene/ncRNA

Inputs

File format	Information contained in file	File description	Notes
bam	alignments	Produced by mapping reads to the transcriptome.	Reads are trimmed using a proprietary version of cutAdapt. We map to transcriptome for a better sensitivity (see details in protocol and example).

Outputs

File format	Information contained in file	File description	Notes
bigWig	signal	Normalized RNA-seq signal	Signals are generated for transcriptome both the plus and minus strands and for unique reads and unique+multimapping reads.
tsv	gene (ncRNA) quantifications	Non-normalized counts.

Running Scripts

Software/Tools

• RSEM

• homer

• HTseq

• FeatureCounts

Example of single case

Set Parameters

# genome seuqneces and annotaions
ln -s /BioII/lulab_b/shared/genomes ~/genomes
export hg38=~/genomes/human_hg38/sequence/GRCh38.p12.genome.fa
export gtf=~/genomes/human_hg38/gtf

# working space
cd ~/proj_exRNA

0.Convert & View

a. homer

b. HTseq

c.FeatureCounts

Example of batch job

Tips/Utilities

Merge multiple bams to one

# mutiple bams to one bam file, in Snakemake form

rule mergeBam:
    input:
        bam = "02.mapping/{library}/sequentialMap/{library}.hg38other.bam"
    output:
        bam = "02.mapping/{library}/sequentialMap/{library}.sequentialMap.bam"
    params:
        jobname = "{library}.mergeBam",
        dir = "02.mapping/{library}/sequentialMap/"
    shell:
        """
        samtools merge {output.bam} {input.bam} `find {params.dir} -name "*.rsem.clean.bam"`

Homework and more

1. Visualize your mapped reads with IGV (locally) and/or UCSC Genome Browser (on line).

2. Learn how to construct the expression matrix using HTSeq, featureCounts and homer; then compare the difference among these three methods.

More Reading and Practice

• Additional Tutorial : 2. Construction of expression matrix

• Bioinformatics Data Skills

  • 11) Working with Alignment Data

Video

a) Expression matrix

@Youtube

@Bilibili