RNA-Seq: A Probabilistic Framework for Aligning Paired-end RNA-seq Data.

[Newsbot!]

Syndicated from PubMed RSS Feeds

A Probabilistic Framework for Aligning Paired-end RNA-seq Data.

Bioinformatics. 2010 Jun 23;

Authors: Hu Y, Wang K, He X, Chiang DY, Prins JF, Liu J

MOTIVATION: The RNA-seq paired-end read (PER) protocol samples transcript fragments longer than the sequencing capability of today's technology by sequencing just the two ends of each fragment. Deep sampling of the transcriptome using the PER protocol presents the opportunity to reconstruct the unsequenced portion of each transcript fragment using end reads from overlapping PERs, guided by the expected length of the fragment. METHODS: A probabilistic framework is described to predict the alignment to the genome of all PER transcript fragments in a PER dataset. Starting from possible exonic and spliced alignments of all end reads, our method constructs potential splicing paths connecting paired ends. An expectation maximization method assigns likelihood values to all splice junctions and assigns the most probable alignment for each transcript fragment. RESULTS: The method was applied to 2x35bp PER datasets from cancer cell lines MCF-7 and SUM-102. PER fragment alignment increased the coverage three fold compared to the alignment of the end reads alone, and increased the accuracy of splice detection. The accuracy of the EM algorithm in the presence of alternative paths in the splice graph was validated by qRT-PCR experiments on 8 exon skipping alternative splicing events. PER fragment alignment with long range splicing confirmed 8 out of 10 fusion events identified in the MCF-7 cell line in an earlier study by Maher et al. (2009). AVAILABILITY: Software available at http://www.netlab.uky.edu/p/bioinfo/PERalign CONTACT: [email protected].

PMID: 20576625 [PubMed - as supplied by publisher]



More...

[repinementer]

nothing found

the link you provide and the author himself provided are not working.

[AndrewCarpenter]

http://bioinformatics.oxfordjournals.org/cgi/content/full/btq336/DC1

[repinementer]

not the journal, SOFTWARE-PERAlign?

not the journal, SOFTWARE-PERAlign?

I already read the articles. what I meant was there is nothing in the link they provided.

[brentp]

this link works, i guessed it from teh parent page:

http://www.netlab.uky.edu/p/bioinfo/MapPER

[repinementer]

MapSplicer Error

When I tried to ran MapSplicer I got this error. COuld any one assist me please. Thanx.

MapSplice_1.13.6]$ python bin/mapsplice_segments.py MapSplice.cfg 

[Fri Aug  6 15:30:44 2010] Beginning Mapsplice run (v1.13.6)
-----------------------------------------------
bin directory: [/home/bogugk/XXX/MapSplice_1.13.6/bin/] 
[Fri Aug  6 15:30:44 2010] Preparing output location /home/bogugk/XXX/MapSplice_1.13.6/
[Fri Aug  6 15:30:44 2010] Checking for chromosomes files or directory
[Fri Aug  6 15:30:44 2010] Checking for Bowtie index files
[Fri Aug  6 15:30:44 2010] Indexing known splices
Traceback (most recent call last):
  File "bin/mapsplice_segments.py", line 5823, in ?
    sys.exit(main())
  File "bin/mapsplice_segments.py", line 5138, in main
    check_bowtie_index(bwt_idx_prefix, chrom_dir_files, 0, FASTA_file_extension)
  File "bin/mapsplice_segments.py", line 3212, in check_bowtie_index
    idx_prefix = build_juncs_bwt_index(chromo_sequences, idx_prefix);
  File "bin/mapsplice_segments.py", line 3189, in build_juncs_bwt_index
    exit(1)
TypeError: 'str' object is not callable

Here is my cfg file

MapSplice_1.13.6]$ cat MapSplice.cfg 
######################################################################
#Configration file for running MapSplice
#Configrate file can be used individually to run MapSplice: python mapsplice_segments.py config_file
#Or mixed with inputs and options: python mapsplice_segments.py [inputs|options] config_file
#Lines start with # are comments
#Some options value can be 'yes' or 'no'
#Values and name are separated by ' = '


######################################################################
#Inputs and outputs
#

#
#Comma separated list of  the RNA-seq read files
#For paired-end reads, the order should be in reads1_end1,reads1_end2,reads2_end1,read2_end2...

reads_file = /home/bogugk/rnaseq/Human_CD4TCell_RNASEQ_data/SRR017234_1.fastq 


#
#The directory containing sequences files of reference genome in FASTA format
#One chromosome per file
#The chromosome name after '>' should not contain tab or blank space
#The chromosome name should be the same of basename of chromosome file
#Suffix of chromosome file name should be 'fa'
#For example if chromosome name after '>' is 'chr1', then the file name should be 'chr1.fa'

chromosome_files_dir = ../MapSplice_1.13.6/bin/bowtie-0.12.5/indexes


#
#The path and base name of index to be searched by Bowtie. 
#If the index does not exist, it will be built from reference genomes indicated by option -c with bowtie-build. 
#The built index will be at bowtie_index_path/bowtie_index_base

Bowtieidx = /home/bogugk/XXX/MapSplice_1.13.6/bin/bowtie-build/
 

#
#The name of the directory in which MapSplice will write its output 

output_dir = /home/bogugk/XXX/MapSplice_1.13.6

#
#Not search alignments in the specified regions in gff format
#(optional)

#avoid_regions = avoid_regions.gff


#
#Search alignments in the specified regions in gff format
#(optional)

#interested_regions = interested_regions.gff


######################################################################
#Basic options
#

#
#Format of input reads, FASTA OR FASTQ
#Reads name after '>' or '@' should not contain blank space or tab

reads_format = FASTQ


#
#If input reads are paired end reads or single reads

#paired_end = yes


#
#Input reads length, read length can be as as short as 36bp, or arbitrary long

read_length = 35


#
#Length of segments reads divided by
#Should be in range [18,25], and should not be longer than half of read length

segment_length = 18


######################################################################
#Advanced options
#

#
#If output semi-canonical and non-canonical junction
#non-canonical: output all junctions include non-canonical, semi-canonical, canonical junctions
#semi-canonical: output semi-canonical, canonical junctions
#canonical:output canonical junctions

junction_type = canonical


#
#If run remap step to increase junction coverage

full_running = yes


#
#The anchor length used for single anchored spliced alignment
#Decrease it will find more alignments but use more running time

anchor_length = 8


#
#If delete temp files to save disk space after finished

remove_temp_files = yes


#
#Number of mismatches allowed in segment mapping
#Should be in range [0,3]

segment_mismatches = 1

#
#Number of mismatches allowed in missed segment(spliced segment)
splice_mismatches = 1


#
#Number of mismatches allowed in remapping step

remap_mismatches = 2


#
#Minimal intron length, the default is 10

min_intron_length = 10


#
#Maximal intron length, the default is 200,000

max_intron_length = 200000


#
#Number of threads for segment mapping, will be used for MapSplice in future

threads = 1


#
#If search whole chromosomes to find splice alignment, instead of searching exonic regions.
#Can find splice alignments in small exons(< segment length), but will use more running time

search_whole_chromosome = no


#
#If yes, try to find spliced alignments and unspliced alignments of a read, and then select best alignment. Will use more running time
#If not, try to find unspliced alignments of a read, if no unspliced alignments found, then try to find spliced alignment

map_segment_directly = no


#
#Run MapPER and generate reads mapping based on probabilistic framework
#More information about probabilistic framework is at bioinformatics.oxfordjournals.org/cgi/reprint/btq336v1.pdf

run_MapPER = no


#
#If output fusion junction, reads should be long enough to be divide more than 2 segments for fusion alignment
#The output is "fusion.junction" and "fusion_junction.unique"

do_fusion = no


#
#Use paired-end reads to generate cluster regions for fusion reads mapping
#Only valid for paired-end reads and full running model(set full_running = yes)
#The output is "cluster.fusion.mapped" and "cluster.fusion.junction"

#do_cluster = no

[raela]

I don't know if this is relevant, but do you already have a built bowtie index?
Bowtieidx = /home/bogugk/XXX/MapSplice_1.13.6/bin/bowtie-build/
^ they want the path to the index, not to bowtie-build, as in something like
../MapSplice_1.13.6/bin/bowtie-0.12.5/indexes/ewbt_base

[greenshell]

Bowtieidx should be the path and base name, not only path
Bowtieidx = /home/bogugk/XXX/MapSplice_1.13.6/bin/bowtie-build/

bowtie and bowtie-build 0.10.1 should be in MapSplice/bin directory

When I tried to ran MapSplicer I got this error. COuld any one assist me please. Thanx.

MapSplice_1.13.6]$ python bin/mapsplice_segments.py MapSplice.cfg 

[Fri Aug  6 15:30:44 2010] Beginning Mapsplice run (v1.13.6)
-----------------------------------------------
bin directory: [/home/bogugk/XXX/MapSplice_1.13.6/bin/] 
[Fri Aug  6 15:30:44 2010] Preparing output location /home/bogugk/XXX/MapSplice_1.13.6/
[Fri Aug  6 15:30:44 2010] Checking for chromosomes files or directory
[Fri Aug  6 15:30:44 2010] Checking for Bowtie index files
[Fri Aug  6 15:30:44 2010] Indexing known splices
Traceback (most recent call last):
  File "bin/mapsplice_segments.py", line 5823, in ?
    sys.exit(main())
  File "bin/mapsplice_segments.py", line 5138, in main
    check_bowtie_index(bwt_idx_prefix, chrom_dir_files, 0, FASTA_file_extension)
  File "bin/mapsplice_segments.py", line 3212, in check_bowtie_index
    idx_prefix = build_juncs_bwt_index(chromo_sequences, idx_prefix);
  File "bin/mapsplice_segments.py", line 3189, in build_juncs_bwt_index
    exit(1)
TypeError: 'str' object is not callable

Here is my cfg file

MapSplice_1.13.6]$ cat MapSplice.cfg 
######################################################################
#Configration file for running MapSplice
#Configrate file can be used individually to run MapSplice: python mapsplice_segments.py config_file
#Or mixed with inputs and options: python mapsplice_segments.py [inputs|options] config_file
#Lines start with # are comments
#Some options value can be 'yes' or 'no'
#Values and name are separated by ' = '


######################################################################
#Inputs and outputs
#

#
#Comma separated list of  the RNA-seq read files
#For paired-end reads, the order should be in reads1_end1,reads1_end2,reads2_end1,read2_end2...

reads_file = /home/bogugk/rnaseq/Human_CD4TCell_RNASEQ_data/SRR017234_1.fastq 


#
#The directory containing sequences files of reference genome in FASTA format
#One chromosome per file
#The chromosome name after '>' should not contain tab or blank space
#The chromosome name should be the same of basename of chromosome file
#Suffix of chromosome file name should be 'fa'
#For example if chromosome name after '>' is 'chr1', then the file name should be 'chr1.fa'

chromosome_files_dir = ../MapSplice_1.13.6/bin/bowtie-0.12.5/indexes


#
#The path and base name of index to be searched by Bowtie. 
#If the index does not exist, it will be built from reference genomes indicated by option -c with bowtie-build. 
#The built index will be at bowtie_index_path/bowtie_index_base

Bowtieidx = /home/bogugk/XXX/MapSplice_1.13.6/bin/bowtie-build/
 

#
#The name of the directory in which MapSplice will write its output 

output_dir = /home/bogugk/XXX/MapSplice_1.13.6

#
#Not search alignments in the specified regions in gff format
#(optional)

#avoid_regions = avoid_regions.gff


#
#Search alignments in the specified regions in gff format
#(optional)

#interested_regions = interested_regions.gff


######################################################################
#Basic options
#

#
#Format of input reads, FASTA OR FASTQ
#Reads name after '>' or '@' should not contain blank space or tab

reads_format = FASTQ


#
#If input reads are paired end reads or single reads

#paired_end = yes


#
#Input reads length, read length can be as as short as 36bp, or arbitrary long

read_length = 35


#
#Length of segments reads divided by
#Should be in range [18,25], and should not be longer than half of read length

segment_length = 18


######################################################################
#Advanced options
#

#
#If output semi-canonical and non-canonical junction
#non-canonical: output all junctions include non-canonical, semi-canonical, canonical junctions
#semi-canonical: output semi-canonical, canonical junctions
#canonical:output canonical junctions

junction_type = canonical


#
#If run remap step to increase junction coverage

full_running = yes


#
#The anchor length used for single anchored spliced alignment
#Decrease it will find more alignments but use more running time

anchor_length = 8


#
#If delete temp files to save disk space after finished

remove_temp_files = yes


#
#Number of mismatches allowed in segment mapping
#Should be in range [0,3]

segment_mismatches = 1

#
#Number of mismatches allowed in missed segment(spliced segment)
splice_mismatches = 1


#
#Number of mismatches allowed in remapping step

remap_mismatches = 2


#
#Minimal intron length, the default is 10

min_intron_length = 10


#
#Maximal intron length, the default is 200,000

max_intron_length = 200000


#
#Number of threads for segment mapping, will be used for MapSplice in future

threads = 1


#
#If search whole chromosomes to find splice alignment, instead of searching exonic regions.
#Can find splice alignments in small exons(< segment length), but will use more running time

search_whole_chromosome = no


#
#If yes, try to find spliced alignments and unspliced alignments of a read, and then select best alignment. Will use more running time
#If not, try to find unspliced alignments of a read, if no unspliced alignments found, then try to find spliced alignment

map_segment_directly = no


#
#Run MapPER and generate reads mapping based on probabilistic framework
#More information about probabilistic framework is at bioinformatics.oxfordjournals.org/cgi/reprint/btq336v1.pdf

run_MapPER = no


#
#If output fusion junction, reads should be long enough to be divide more than 2 segments for fusion alignment
#The output is "fusion.junction" and "fusion_junction.unique"

do_fusion = no


#
#Use paired-end reads to generate cluster regions for fusion reads mapping
#Only valid for paired-end reads and full running model(set full_running = yes)
#The output is "cluster.fusion.mapped" and "cluster.fusion.junction"

#do_cluster = no

[repinementer]

yes. It ran well after correcting it but only up to certain extent. I gave up. TopHat is working fine for me now. But thanx for the post above

[greenshell]

The segment length should not exceed half of read length, try set segment_length = 17 for read length 35

When I tried to ran MapSplicer I got this error. COuld any one assist me please. Thanx.

MapSplice_1.13.6]$ python bin/mapsplice_segments.py MapSplice.cfg 

[Fri Aug  6 15:30:44 2010] Beginning Mapsplice run (v1.13.6)
-----------------------------------------------
bin directory: [/home/bogugk/XXX/MapSplice_1.13.6/bin/] 
[Fri Aug  6 15:30:44 2010] Preparing output location /home/bogugk/XXX/MapSplice_1.13.6/
[Fri Aug  6 15:30:44 2010] Checking for chromosomes files or directory
[Fri Aug  6 15:30:44 2010] Checking for Bowtie index files
[Fri Aug  6 15:30:44 2010] Indexing known splices
Traceback (most recent call last):
  File "bin/mapsplice_segments.py", line 5823, in ?
    sys.exit(main())
  File "bin/mapsplice_segments.py", line 5138, in main
    check_bowtie_index(bwt_idx_prefix, chrom_dir_files, 0, FASTA_file_extension)
  File "bin/mapsplice_segments.py", line 3212, in check_bowtie_index
    idx_prefix = build_juncs_bwt_index(chromo_sequences, idx_prefix);
  File "bin/mapsplice_segments.py", line 3189, in build_juncs_bwt_index
    exit(1)
TypeError: 'str' object is not callable

Here is my cfg file

MapSplice_1.13.6]$ cat MapSplice.cfg 
######################################################################
#Configration file for running MapSplice
#Configrate file can be used individually to run MapSplice: python mapsplice_segments.py config_file
#Or mixed with inputs and options: python mapsplice_segments.py [inputs|options] config_file
#Lines start with # are comments
#Some options value can be 'yes' or 'no'
#Values and name are separated by ' = '


######################################################################
#Inputs and outputs
#

#
#Comma separated list of  the RNA-seq read files
#For paired-end reads, the order should be in reads1_end1,reads1_end2,reads2_end1,read2_end2...

reads_file = /home/bogugk/rnaseq/Human_CD4TCell_RNASEQ_data/SRR017234_1.fastq 


#
#The directory containing sequences files of reference genome in FASTA format
#One chromosome per file
#The chromosome name after '>' should not contain tab or blank space
#The chromosome name should be the same of basename of chromosome file
#Suffix of chromosome file name should be 'fa'
#For example if chromosome name after '>' is 'chr1', then the file name should be 'chr1.fa'

chromosome_files_dir = ../MapSplice_1.13.6/bin/bowtie-0.12.5/indexes


#
#The path and base name of index to be searched by Bowtie. 
#If the index does not exist, it will be built from reference genomes indicated by option -c with bowtie-build. 
#The built index will be at bowtie_index_path/bowtie_index_base

Bowtieidx = /home/bogugk/XXX/MapSplice_1.13.6/bin/bowtie-build/
 

#
#The name of the directory in which MapSplice will write its output 

output_dir = /home/bogugk/XXX/MapSplice_1.13.6

#
#Not search alignments in the specified regions in gff format
#(optional)

#avoid_regions = avoid_regions.gff


#
#Search alignments in the specified regions in gff format
#(optional)

#interested_regions = interested_regions.gff


######################################################################
#Basic options
#

#
#Format of input reads, FASTA OR FASTQ
#Reads name after '>' or '@' should not contain blank space or tab

reads_format = FASTQ


#
#If input reads are paired end reads or single reads

#paired_end = yes


#
#Input reads length, read length can be as as short as 36bp, or arbitrary long

read_length = 35


#
#Length of segments reads divided by
#Should be in range [18,25], and should not be longer than half of read length

segment_length = 18


######################################################################
#Advanced options
#

#
#If output semi-canonical and non-canonical junction
#non-canonical: output all junctions include non-canonical, semi-canonical, canonical junctions
#semi-canonical: output semi-canonical, canonical junctions
#canonical:output canonical junctions

junction_type = canonical


#
#If run remap step to increase junction coverage

full_running = yes


#
#The anchor length used for single anchored spliced alignment
#Decrease it will find more alignments but use more running time

anchor_length = 8


#
#If delete temp files to save disk space after finished

remove_temp_files = yes


#
#Number of mismatches allowed in segment mapping
#Should be in range [0,3]

segment_mismatches = 1

#
#Number of mismatches allowed in missed segment(spliced segment)
splice_mismatches = 1


#
#Number of mismatches allowed in remapping step

remap_mismatches = 2


#
#Minimal intron length, the default is 10

min_intron_length = 10


#
#Maximal intron length, the default is 200,000

max_intron_length = 200000


#
#Number of threads for segment mapping, will be used for MapSplice in future

threads = 1


#
#If search whole chromosomes to find splice alignment, instead of searching exonic regions.
#Can find splice alignments in small exons(< segment length), but will use more running time

search_whole_chromosome = no


#
#If yes, try to find spliced alignments and unspliced alignments of a read, and then select best alignment. Will use more running time
#If not, try to find unspliced alignments of a read, if no unspliced alignments found, then try to find spliced alignment

map_segment_directly = no


#
#Run MapPER and generate reads mapping based on probabilistic framework
#More information about probabilistic framework is at bioinformatics.oxfordjournals.org/cgi/reprint/btq336v1.pdf

run_MapPER = no


#
#If output fusion junction, reads should be long enough to be divide more than 2 segments for fusion alignment
#The output is "fusion.junction" and "fusion_junction.unique"

do_fusion = no


#
#Use paired-end reads to generate cluster regions for fusion reads mapping
#Only valid for paired-end reads and full running model(set full_running = yes)
#The output is "cluster.fusion.mapped" and "cluster.fusion.junction"

#do_cluster = no

[rocksd]

Hi greenshell,

I am about to use MapSplice as I want to look at non-canonical splice sites. But I have a few questions:

1. Do I have to build Bowtie index for every chr individually? I guess the answer is yes but just want to confirm. I have a Bowtie index that I built for Tophat, but it is for the whole genome.

2. I am a little confused with your definition of paired-end reads and ends. If I have 8 different samples for a paired-end run, each occupying one lane in GAIIX, so I will have a total 16 read files, which will be named like this: s_1_reads1, s_1_reads2, s_2_read1, s_2_reads2, s_3_reads1, s_3_reads2, etc where "s_1" means "sample run on lane 1". Since I have 8 different samples, I will have to map one sample at a time. In this case, should I use "reads1_end1,reads1_end2" or "read_base_name/1 and read_base_name/2" for -u/ option?

Thanks in advance!

[greenshell]

1. Do I have to build Bowtie index for every chr individually? I guess the answer is yes but just want to confirm. I have a Bowtie index that I built for Tophat, but it is for the whole genome.

No, the bowtie index should be built for whole chromosome. But the chromosome files should contain 1 chromosome per file

2. I am a little confused with your definition of paired-end reads and ends. If I have 8 different samples for a paired-end run, each occupying one lane in GAIIX, so I will have a total 16 read files, which will be named like this: s_1_reads1, s_1_reads2, s_2_read1, s_2_reads2, s_3_reads1, s_3_reads2, etc where "s_1" means "sample run on lane 1". Since I have 8 different samples, I will have to map one sample at a time. In this case, should I use "reads1_end1,reads1_end2" or "read_base_name/1 and read_base_name/2" for -u/ option?

You should use "reads1_end1,reads1_end2" for -u option. You can run the 16 read files together in the way "reads1_end1,reads1_end2,reads2_end1,reads2_end2.....reads8_end1,reads8_end2"

In the two read ends files.
If the first read is
@HWI-EAS7_0001:1:1:3:1132#0/1
in one read end file.

The the first read in the other read end file should be
@HWI-EAS7_0001:1:1:3:1132#0/2