Skip to content

Latest commit

 

History

History
304 lines (230 loc) · 12 KB

1.mapping-annotation-and-qc.md

File metadata and controls

304 lines (230 loc) · 12 KB
Raw

1.Mapping, Annotation and QC

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

Inputs

File format Information contained in file File description Notes
fastq reads small RNA-seq reads from single-end libraries. This will be processed from long reads we got. Raw reads are trimmed as single-ended, 50nt clean reads (see details in the protocol and example).
fasta genome/transcriptome sequences sequences of the whole human genome or different RNAs (transcriptome)
gtf / gff genome annotation Default genome annotation file is from GENCODE, mirBase, etc. Be careful about exon and intron for mRNA/lncRNA, precursor and mature product for miRNA.

Outputs

| File format | Information contained in file | File description | Notes | | --- | --- | --- | | sam/bam | alignments, i.e. mapped reads (bam is a binary format of sam) | Produced by mapping reads to transcriptome | We map to transcriptome for a better sensitivity (see details in protocol and example). | | bt2 | indexed sequences | Indexed genome/ transcriptome | Different ncRNA types are indexed separately and mapped sequentially. |

Running Scripts

Software/Tools

  • RSEM: index and many utility package for transcriptome analysis
  • fastqc: QC for fastq file
  • cutadapt: trim and cut adaptor
  • bowtie2: mapping

Example of single case

{% tabs %} {% tab title="Set Parameters" %}

# genome seuqneces and annotaions
ln -s /BioII/lulab_b/shared/genomes $HOME/genomes
export hg38=$HOME/genomes/human_hg38/sequence/GRCh38.p10.genome.fa
export gtf=$HOME/genomes/human_hg38/anno/gtf
export index=$HOME/genomes/human_hg38/index/transcriptome_rsem_bowtie2

# raw data (fastq files)
cd $HOME/proj_exRNA/example/fastq
ln -s /BioII/lulab_b/shared/projects/exRNA/example/100K/*.fastq .

# working space
cd $HOME/proj_exRNA

{% endtab %}

{% tab title="0.Index" %}

# Produce indexed transcriptome, using miRNA as an example
# Usually, you do not need to do this because we have already index one
# in "$index"
# We need RSEM to index transcriptome (RNAs) because we'll convert 
# the relative coordinates to genomic coordinats after mapping to RNAs

#############################################################
# 0: build bowtie2 index using RSEM for RNA transcriptome
#############################################################

rsem-prepare-reference --gtf $gtf/miRNA.gencode27.gtf --bowtie2 $hg38 \
  RNA_index/miRNA.gencode27 > log.0 2>err.0 &
  # RSEM will take exons only 
  # only if you prepare gtf files in the right format

{% endtab %}

{% tab title="1.1QC-1.2 Trim-1.3QC" %}

# check raw reads' quality
# trim low quality ends + keep first 50nt + remove 3' adapter

#############################################################
# 1.1: fastaqc - repeat fastqc before and after each trim step
#############################################################

fastqc -q -o example/fastqc NC_1.fastq \
  > log.1.1 2>err.1.1 &
    #output the QC files to a dir examples/fastqc

#############################################################
# 1.2: trim - cut adaptor + trim long read
#############################################################
cutadapt -u -100 -q 30,30 --trim-n -m 15 -a \
 AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -o NC_1.fastq.trim NC_1.fastq  
    # -u -100    remove last 100nt so that the first 50nt is kept
    # -q 30,30   read quality need to be above 30
    # -m 15      reads less than 15nt are removed


#############################################################
# repeat 1.3 QC after 1.2 Trim
# make sure the low quality reads have been removed and/or trimmed.
#############################################################

fastqc -q -o example/fastqc NC_1.fastq.trim \
  > log.1.3 2>err.1.3 &

{% endtab %}

{% tab title="1.4 Clean rRNA reads" %}

############################################################# 
# 1.4: clean rRNA reads 
#############################################################

cd $proj/fastq 
for i in `ls *.fastq` ; do
	sample="${i/.fastq/}"
	sam=$sample.rRNA.sam
	bam=$sample.rRNA.bam
	fqin=$sample.fastq.trim
	fqout=$sample.no_rRNA.fq
	RNA=$index/rRNA
	echo "start $sample :"
		
	bowtie2 -p 4 --sensitive-local --norc --no-unal --un ../mapped/$fqout \
     -x $RNA $fqin -S ../mapped/$sam 
	#map to rRNAs (exons only, sense strand only)
	# (default) look for multiple alignments, report best, with MAPQ
	# --sensitive-local(default): allow no mismatch, etc
	# --norc:  do not align reverse-complement version of read
	# --no-unal: suppress SAM records for unaligned reads
	# --un: store unmapped reads
	# -x: indexed genome/transcriptome 
	# -S: output file foramt as sam

	samtools view -b ../mapped/$sam > ../mapped/$bam
	#convert sam to bam for mapped reads on rRNAs (sense strand)

	echo "$sample finished."
done

{% endtab %}

{% tab title="1.5 Mapping" %}

# map to different RNAs (transcriptome) step by step
#############################################################
# 1.5 align cleaned reads to transcriptome (different RNAs) 
# first align to miRNAs
#############################################################
cd $proj/fastq 
for i in `ls *.fastq`  ; do
	sample="${i/.fastq/}"
	fqin=$sample.no_rRNA.fq
	fqout=$sample.no_miRNA.fq
	sam=$sample.miRNA.sam
	bam=$sample.miRNA.bam
	RNA=$index/miRNA
	echo "start $sample :"
	
	bowtie2 -p 4 --sensitive-local --norc --no-unal --un ../mapped/$fqout \
    	 -x $RNA ../mapped/$fqin -S ../mapped/$sam 
	#map to miRNAs (pre-miRNA or pri-miRNA, sense strand only)
	# (default) look for multiple alignments, report best, with MAPQ
	# --sensitive-local(default): allow no mismatch, etc
	# --norc:  do not align reverse-complement version of read
	# --no-unal: suppress SAM records for unaligned reads
	# --un: store unmapped reads
	# -x: indexed genome/transcriptome 
	# -S: output file foramt as sam

	samtools view -b ../mapped/$sam > ../mapped/$bam
	#convert sam to bam for mapped reads on miRNAs (sense strand)
	
	echo "$sample finished."
done


##############################################################
# then align to next RNA (e.g. piRNA)
# repeat mapping ...

{% endtab %} {% endtabs %}

Example of batch job

Tips/Utilities

A better view of gff/gtf file

less -S *.gff  # chop long lines rather than wrap them

Convert gff to gtf

gffread tRNA.gff -T -o tRNA.gtf

Lift over gff/bed/bam version

# convert bed
liftOver old.bed map.chain new.bed unmapped

# convert gtf
## option 1: simple method (not recommended)
liftOver -gff old.gtf map.chain new.gtf unmapped
###   -gff  File is in gff/gtf format.  Note that the gff lines are converted
###         separately.  It would be good to have a separate check after this
###         that the lines that make up a gene model still make a plausible gene
###         after liftOver

## option 2: conver to genePred format first
gtfToGenePred -genePredExt old.gtf old.gp
liftOver -genePred old.gp map.chain new.gp unmapped
genePredToGtf file new.gp new.gtf

Collapse duplicated reads

Not needed for sRNA-seq.

#option 1: collapse before mapping
fastx_collapser -Q 33 -i input.fastq -o output.collapsed.fastq

#option 2: Collapse after mapping
samtools rmdup -s sorted.bam output.dedup.bam

java -jar picard.jar MarkDuplicates \
I=input.bam \
O=output.markDuplicates \
M=output.metric \
REMOVE_DUPLICATES=true

Homework and more

1 Map on example data

Map the example fastq files and summarize the distribution of 1) reads' length and 2) reads' percentage for each type of RNA (e.g. miRNA, piRNA, Y RNA, srp RNA, etc) ****. (You may find some script examples for this statistics calculation from this link.)

2 Sample QC

Use the following criteria to QC the mapped samples for cfRNAs (exRNAs):

| Check point | Threshold | Notes | | --- | --- | --- | --- | --- | --- | --- | --- | | Raw reads quality | reads quality >28 (median lines in green area) | Check fastqc results(*.html) | | Clean reads number | > 10 million | Adaptors and too-short sequences removed reads | | rRNAs% | < 10% | Reads mapped to rRNAs (all following % are divided by the total number of clean reads) | | HG% | > 60% (optional) | Reads mapped to Human Genome except rRNAs | | Transcriptome% | > 50% | Reads mapped to Human Transcriptome (including rRNA, miRNA, piRNA, Y RNA, srpRNA, snRNA, snoRNA, tRNA, mRNA exons, lncRNA exons, TUCP exons) | | Y RNA% | 10%~65% | Reads mapped to Y RNA | | miRNA% | 10%~65% (up to 80% for exoRNAs) | Reads mapped to miRNA |

3 A fast multi mapping methods to test different mapping orders

3.1 Why use this method: We want to determine a satisfying mapping order by examine the mapping ratio and length distribution. It takes hours to map and collect statistics from a specific order. We use a special method, using python (numpy) to manipulate matrix. By using this method, we can test hundreds of mapping order within few minutes.

3.2 Scripts

At first we do mapping for all RNAs at one time, we use this snakefile and config file to do mapping at one time. (One map snakefile. One map config.)

Secondly we do some matrix manipulation tricks. Using numpy to do matrix indexing and sorting is ultra fast. By using some tricks, we can have the desired statistics within few minutes for hundreds of different mapping orders: (Mapping order for ratio. Mapping order for length distribution)(You can find some tries and tests of the mapping order methods and statistics and plot jupyter notebook here: Statistics and plot)

At last we also extend our algorithm, we can create a big dictionary to store all the possible mapping assignment in every order and every one map condition. We encode the binary condition sequence into 10-bit digits. And store the whole dictionary together. It can test up to 100,000 kinds of mapping order within few hours. The jupyter-notebook presents the main idea of this method. (Ultra-fast mapping methods)

More Reading and Practice