Matt:LabNotes/2014-7-30

From ZhangLabWiki
Jump to navigation Jump to search

RNA-SeQC[edit]

  • Tophat aligned to hg19 BAM input

Trial 1[edit]

Add Read Groups and coordinate sort[edit]

  • Picard.AddOrReplaceReadGroups
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/AddOrReplaceReadGroups.jar INPUT='accepted_hits.bam' OUTPUT='accepted_hits_RanHex.coorsortRG.bam' SORT_ORDER=coordinate RGID='1' RGLB='RTprimer' RGPL='ILLUMINA' RGPU='flowcell-barcode.lane' RGSM='RanHex'

Index reference fasta and create sequence dictionary[edit]

  • reference fasta index already exists
 /home/kunzhang/HsGenome/hg19/HsGenome19.fa.fai
  • Picard.CreateSequenceDictionary
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/CreateSequenceDictionary.jar REFERENCE=/home/kunzhang/HsGenome/hg19/HsGenome19.fa OUTPUT=/home/mzcai/Genomes/HsGenome19.dict

Match order of contigs in Bam file to reference file[edit]

  • Picard.ReorderSam
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/ReorderSam.jar INPUT=accepted_hits_RanHex.coorsortRG.bam OUTPUT=accepted_hits_RanHex.coorsortRGreorder.bam REFERENCE=/home/mzcai/Genomes/HsGenome19.fa

Index bam file[edit]

 samtools index accepted_hits_RanHex.coorsortRGreorder.bam

Create sample file[edit]

  • ~/scratch/RanHex_EnrichmRNA/inVitroRTSeq/RNA-SeQC_samples.txt
 Sample ID       Bam File        Notes
 RanHex  /home/mzcai/scratch/RanHex_EnrichmRNA/inVitroRTSeq/tophat_hg19unmask_Indx26_RanHex/accepted_hits_RanHex.coorsortRGreorder.bam   Indx26RanHex

Run RNA-SeQC v1.1.7.jar[edit]

 java -jar /home/kunzhang/softwares/RNA-SeQC_v1.1.7.jar -bwa /home/kunzhang/softwares/bwa-0.7.5a/bwa -BWArRNA /home/mzcai/scratch/RanHex_EnrichmRNA/mart_export_exons_totalrRNA.fa -o /home/mzcai/scratch/RanHex_EnrichmRNA/inVitroRTSeq/RNA-SeQC_output -r /home/mzcai/Genomes/HsGenome19.fa -s /home/mzcai/scratch/RanHex_EnrichmRNA/inVitroRTSeq/RNA-SeQC_samples.txt -singleEnd -t /home/mzcai/LTS/Genomes/Homo_sapiens.GRCh37.75.chr.gtf -ttype 2
 The required transcript_id attribute was not found on line chr1 pseudogene      gene    11869   14412   .       +       .       gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene";

Trial 2[edit]

  • Use provided gencode gtf annotation file
    • Also use provided rRNA ref file and GC definitions file to be consistent with example

Download resources from broadinstitute[edit]

  • ~/Genomes/RNA-SeQCResources
  • Download human rRNA reference file
 human_all_rRNA.fasta
  • Download GENCODE GTF file
 gencode.v7.annotation.gtf
  • Download GENCODE GC definitions file
 gencode.v7.gc.txt

Run RNA-SeQC v1.1.7.jar[edit]

 java -jar /home/kunzhang/softwares/RNA-SeQC_v1.1.7.jar -bwa /home/kunzhang/softwares/bwa-0.7.5a/ -BWArRNA ~/Genomes/RNA-SeqCResources/human_all_rRNA.fasta -s ~/scratch/RanHex_EnrichmRNA/inVitroRTSeq/RNA-SeQC_samples.txt -t ~/Genomes/RNA-SeqCResources/gencode.v7.annotation.gtf -r ~/Genomes/HsGenome19.fa -o ~/scratch/RanHex_EnrichmRNA/inVitroRTSeq/RNA-SeQC_output/ -strat gc -gc ~/Genomes/RNA-SeqCResources/gencode.v7.gc.txt -singleEnd
 
 org.broadinstitute.sting.utils.exceptions.UserException$LexicographicallySortedSequenceDictionary: Lexicographically sorted human genome sequence detected in reads.
 For safety's sake the GATK requires human contigs in karyotypic order: 1, 2, ..., 10, 11, ..., 20, 21, 22, X, Y with M either leading or trailing these contigs.
 This is because all distributed GATK resources are sorted in karyotypic order, and your processing will fail when you need to use these files.
 You can use the ReorderSam utility to fix this problem: http://www.broadinstitute.org/gsa/wiki/index.php/ReorderSam

Trial 3[edit]

Reorder HsGenome19.fa to karyotypic order[edit]

 csplit -f c HsGenome19.fa /'>chr9'/ /'>chr5'/ /'>chr2'/ /'>chr17'/ /'>chr14'/ /'>chr4'/ /'>chr16'/ /'>chr21'/ /'>chrM'/ /'>chr7'/ /'>chr3'/ /'>chr18'/ /'>chr12'/ /'>chrX'/ /'>chr13'/ /'>chr15'/ /'>chr8'/ /'>chr22'/ /'>chr11'/ /'>chr10'/ /'>chr20'/ /'>chr19'/ /'>chr6'/ /'>chr1'/ /'>chrY'/
 cat c24 c03 c11 c06 c02 c23 c10 c17 c01 c20 c19 c13 c15 c05 c16 c07 c04 c12 c22 c21 c08 c18 c14 c25 c09 > HsGenome19_karyorder.fa
 rm HsGenome19.*
 mv HsGenome19_karyorder.fa ./HsGenome19.fa
Old chr order Add // csplit output files cat order New chr order
>chr9 /'>chr9'/ c01 c24 chr1
>chr5 /'>chr5'/ c02 c03 chr2
>chr2 /'>chr2'/ c03 c11 chr3
>chr17 /'>chr17'/ c04 c06 chr4
>chr14 /'>chr14'/ c05 c02 chr5
>chr4 /'>chr4'/ c06 c23 chr6
>chr16 /'>chr16'/ c07 c10 chr7
>chr21 /'>chr21'/ c08 c17 chr8
>chrM /'>chrM'/ c09 c01 chr9
>chr7 /'>chr7'/ c10 c20 chr10
>chr3 /'>chr3'/ c11 c19 chr11
>chr18 /'>chr18'/ c12 c13 chr12
>chr12 /'>chr12'/ c13 c15 chr13
>chrX /'>chrX'/ c14 c05 chr14
>chr13 /'>chr13'/ c15 c16 chr15
>chr15 /'>chr15'/ c16 c07 chr16
>chr8 /'>chr8'/ c17 c04 chr17
>chr22 /'>chr22'/ c18 c12 chr18
>chr11 /'>chr11'/ c19 c22 chr19
>chr10 /'>chr10'/ c20 c21 chr20
>chr20 /'>chr20'/ c21 c08 chr21
>chr19 /'>chr19'/ c22 c18 chr22
>chr6 /'>chr6'/ c23 c14 chrX
>chr1 /'>chr1'/ c24 c25 chrY
>chrY /'>chrY'/ c25 c09 chrM

Index reference fasta and create sequence dictionary[edit]

  • samtools faidx
 samtools faidx ~/Genomes/HsGenome19.fa
  • Picard.CreateSequenceDictionary
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/CreateSequenceDictionary.jar REFERENCE=~/Genomes/HsGenome19.fa OUTPUT=~/Genomes/HsGenome19.dict

Match order of contigs in Bam file to reference file[edit]

  • Picard.ReorderSam
 rm accepted_hits_RanHex.coorsortRGreorder.bam*
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/ReorderSam.jar INPUT=accepted_hits_RanHex.coorsortRG.bam OUTPUT=accepted_hits_RanHex.coorsortRGreorder.bam REFERENCE=~/Genomes/HsGenome19.fa

Index bam file[edit]

 samtools index accepted_hits_RanHex.coorsortRGreorder.bam

Run RNA-SeQC v1.1.7.jar[edit]

 java -jar /home/kunzhang/softwares/RNA-SeQC_v1.1.7.jar -bwa /home/kunzhang/softwares/bwa-0.7.5a/bwa -BWArRNA ~/Genomes/RNA-SeqCResources/human_all_rRNA.fasta -s ~/scratch/RanHex_EnrichmRNA/inVitroRTSeq/RNA-SeQC_samples.txt -t ~/Genomes/RNA-SeqCResources/gencode.v7.annotation.gtf -r ~/Genomes/HsGenome19.fa -o ~/scratch/RanHex_EnrichmRNA/inVitroRTSeq/RNA-SeQC_output/ -strat gc -gc ~/Genomes/RNA-SeqCResources/gencode.v7.gc.txt -singleEnd


Trial 4: Add Other Samples[edit]

  • In respective directories
    • eg. ~/scratch/RanHex_EnrichmRNA/inVitroRTSeq/tophat_hg19unmask_Indx27_dT/

Add Read Groups and coordinate sort[edit]

  • Picard.AddOrReplaceReadGroups
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/AddOrReplaceReadGroups.jar INPUT='accepted_hits.bam' OUTPUT='accepted_hits_dT.coorsortRG.bam' SORT_ORDER=coordinate RGID='2' RGLB='RTprimer' RGPL='ILLUMINA' RGPU='flowcell-barcode.lane' RGSM='dT'
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/AddOrReplaceReadGroups.jar INPUT='accepted_hits.bam' OUTPUT='accepted_hits_FISSEQRT.coorsortRG.bam' SORT_ORDER=coordinate RGID='3' RGLB='RTprimer' RGPL='ILLUMINA' RGPU='flowcell-barcode.lane' RGSM='FISSEQRT'
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/AddOrReplaceReadGroups.jar INPUT='accepted_hits.bam' OUTPUT='accepted_hits_Top48.coorsortRG.bam' SORT_ORDER=coordinate RGID='4' RGLB='RTprimer' RGPL='ILLUMINA' RGPU='flowcell-barcode.lane' RGSM='Top48'
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/AddOrReplaceReadGroups.jar INPUT='accepted_hits.bam' OUTPUT='accepted_hits_gDNA.coorsortRG.bam' SORT_ORDER=coordinate RGID='5' RGLB='RTprimer' RGPL='ILLUMINA' RGPU='flowcell-barcode.lane' RGSM='gDNA'

Match order of contigs in Bam file to reference file[edit]

  • Picard.ReorderSam
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/ReorderSam.jar INPUT=accepted_hits_dT.coorsortRG.bam OUTPUT=accepted_hits_dT.coorsortRGreorder.bam REFERENCE=/home/mzcai/Genomes/HsGenome19.fa
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/ReorderSam.jar INPUT=accepted_hits_FISSEQRT.coorsortRG.bam OUTPUT=accepted_hits_FISSEQRT.coorsortRGreorder.bam REFERENCE=/home/mzcai/Genomes/HsGenome19.fa
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/ReorderSam.jar INPUT=accepted_hits_Top48.coorsortRG.bam OUTPUT=accepted_hits_Top48.coorsortRGreorder.bam REFERENCE=/home/mzcai/Genomes/HsGenome19.fa
 java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.79/ReorderSam.jar INPUT=accepted_hits_gDNA.coorsortRG.bam OUTPUT=accepted_hits_gDNA.coorsortRGreorder.bam REFERENCE=/home/mzcai/Genomes/HsGenome19.fa

Index bam file[edit]

 samtools index accepted_hits_dT.coorsortRGreorder.bam
 samtools index accepted_hits_FISSEQRT.coorsortRGreorder.bam
 samtools index accepted_hits_Top48.coorsortRGreorder.bam
 samtools index accepted_hits_gDNA.coorsortRGreorder.bam

Create sample file[edit]

  • ~/scratch/RanHex_EnrichmRNA/inVitroRTSeq/RNA-SeQC_samples.txt
 Sample ID       Bam File        Notes
 RanHex  /home/mzcai/scratch/RanHex_EnrichmRNA/inVitroRTSeq/tophat_hg19unmask_Indx26_RanHex/accepted_hits_RanHex.coorsortRGreorder.bam   Indx26RanHex
 dT      /home/mzcai/scratch/RanHex_EnrichmRNA/inVitroRTSeq/tophat_hg19unmask_Indx27_dT/accepted_hits_dT.coorsortRGreorder.bam   Indx27dT
 FISSEQRT        /home/mzcai/scratch/RanHex_EnrichmRNA/inVitroRTSeq/tophat_hg19unmask_Indx28_FISSEQRT/accepted_hits_FISSEQRT.coorsortRGreorder.bam       Indx28FISSEQRT
 Top48   /home/mzcai/scratch/RanHex_EnrichmRNA/inVitroRTSeq/tophat_hg19unmask_Indx29_Top48/accepted_hits_Top48.coorsortRGreorder.bam     Indx29Top48
 gDNA    /home/mzcai/scratch/RanHex_EnrichmRNA/inVitroRTSeq/tophat_hg19unmask_Indx32_gDNA/accepted_hits_gDNA.coorsortRGreorder.bam       Indx32gDNA

Run RNA-SeQC v1.1.7.jar[edit]

 java -jar /home/kunzhang/softwares/RNA-SeQC_v1.1.7.jar -bwa /home/kunzhang/softwares/bwa-0.7.5a/bwa -BWArRNA ~/Genomes/RNA-SeqCResources/human_all_rRNA.fasta -s ~/scratch/RanHex_EnrichmRNA/inVitroRTSeq/RNA-SeQC_samples.txt -t ~/Genomes/RNA-SeqCResources/gencode.v7.annotation.gtf -r ~/Genomes/HsGenome19.fa -o ~/scratch/RanHex_EnrichmRNA/inVitroRTSeq/RNA-SeQC_outputALL/ -strat gc -gc ~/Genomes/RNA-SeqCResources/gencode.v7.gc.txt -singleEnd
  • rRNA and rRNA rate values for all samples in index.html are "NA"