Matt:LabNotes/2014-7-30: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Mzcai
mNo edit summary
>Mzcai
mNo edit summary
Line 1: Line 1:
==RNA-SeQC==
==RNA-SeQC==
*Tophat aligned to hg19 BAM input


*Bam input
===Tial 1===
**Add Read Groups and '''coordinate sort'''
====Add Read Groups and '''coordinate sort'''====
***Picard.AddOrReplaceReadGroups
*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'
   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
 
***reference fasta index already exists
====Index reference fasta and create sequence dictionary====
*reference fasta index already exists
   /home/kunzhang/HsGenome/hg19/HsGenome19.fa.fai
   /home/kunzhang/HsGenome/hg19/HsGenome19.fa.fai
***Picard.CreateSequenceDictionary
*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
   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
 
***Picard.ReorderSam
====Match order of contigs in Bam file to reference file====
*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
   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
 
====Index bam file====
   samtools index accepted_hits_RanHex.coorsortRGreorder.bam
   samtools index accepted_hits_RanHex.coorsortRGreorder.bam


*Create sample file: ~/scratch/RanHex_EnrichmRNA/inVitroRTSeq/RNA-SeQC_samples.txt
====Create sample file====
*~/scratch/RanHex_EnrichmRNA/inVitroRTSeq/RNA-SeQC_samples.txt
   Sample ID      Bam File        Notes
   Sample ID      Bam File        Notes
   RanHex  /home/mzcai/scratch/RanHex_EnrichmRNA/inVitroRTSeq/tophat_hg19unmask_Indx26_RanHex/accepted_hits_RanHex.coorsortRGreorder.bam  Indx26RanHex
   RanHex  /home/mzcai/scratch/RanHex_EnrichmRNA/inVitroRTSeq/tophat_hg19unmask_Indx26_RanHex/accepted_hits_RanHex.coorsortRGreorder.bam  Indx26RanHex


====Run RNA-SeQC v1.1.7.jar====
   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
   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";
   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";


===Download resources from [http://www.broadinstitute.org/cancer/cga/rnaseqc_download broadinstitute]===
===Trial 2===
*Use provided gencode gtf annotation file
**Also use provided rRNA ref file and GC definitions file to be consistent with example
 
====Download resources from [http://www.broadinstitute.org/cancer/cga/rnaseqc_download broadinstitute]====
*~/Genomes/RNA-SeQCResources
*~/Genomes/RNA-SeQCResources
*Download human rRNA reference file
*Download human rRNA reference file
Line 32: Line 43:
   gencode.v7.gc.txt
   gencode.v7.gc.txt


====Run RNA-SeQC v1.1.7.jar====
   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
   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.
   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.
  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.
  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
  You can use the ReorderSam utility to fix this problem: http://www.broadinstitute.org/gsa/wiki/index.php/ReorderSam


 
===Trial 3===
===Reorder HsGenome19.fa to karyotypic order===
====Reorder HsGenome19.fa to karyotypic order====
   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'/
   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
   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


{| {{table}}
{| {{table}}
Line 101: Line 115:
| >chrY||/'>chrY'/||c25||c09||chrM
| >chrY||/'>chrY'/||c25||c09||chrM
|}
|}
====

Revision as of 17:58, 1 August 2014

RNA-SeQC

  • Tophat aligned to hg19 BAM input

Tial 1

Add Read Groups and coordinate sort

  • 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

  • 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

  • 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

 samtools index accepted_hits_RanHex.coorsortRGreorder.bam

Create sample file

  • ~/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

 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

  • 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

  • ~/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

 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

Reorder HsGenome19.fa to karyotypic order

 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

==