Query with Biostrings pairwise alignment output positon
1
0
Entering edit mode
dhwani2410 • 0
@dhwani2410-14626
Last seen 6.0 years ago
> query <- readDNAStringSet("ref.fasta")
> sub1 <- readDNAStringSet("seq1.fasta")
> sub2 <- readDNAStringSet("seq2.fasta")
> mat <- nucleotideSubstitutionMatrix(match = 5, mismatch = -4, baseOnly = TRUE)
> res1 <- pairwiseAlignment(pattern = query,subject = sub1,gapOpening = 10.5, gapExtension = .5,substitutionMatrix = mat)
> res2 <- pairwiseAlignment(pattern = query,subject = sub2,gapOpening = 10.5, gapExtension = .5,substitutionMatrix = mat)
> deletion(res1)
IRangesList of length 1
[[1]]
IRanges of length 1
    start end width
[1]   533 533     1

> insertion(res2);deletion(res2)
IRangesList of length 1
[[1]]
IRanges of length 1
    start end width
[1]   363 364     2

IRangesList of length 1
[[1]]
IRanges of length 3
    start end width
[1]   196 213    18
[2]   201 203     3
[3]   336 337     2

> writePairwiseAlignments(res1)
########################################
# Program: Biostrings (version 2.38.4), a Bioconductor package
# Rundate: Thu Dec 14 11:35:27 2017
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: ref
# 2: seq1
# Matrix: NA
# Gap_penalty: 11.0
# Extend_penalty: 0.5
#
# Length: 1099
# Identity:    1098/1099 (99.9%)
# Similarity:    NA/1099 (NA%)
# Gaps:           1/1099 (0.1%)
# Score: 5479
#
#
#=======================================

ref                1 ATGGCCGTCATGGCGCCCCGAACCCTCCTCCTGCTACTCTCGGGGGCCCT     50
                     ||||||||||||||||||||||||||||||||||||||||||||||||||
seq1               1 ATGGCCGTCATGGCGCCCCGAACCCTCCTCCTGCTACTCTCGGGGGCCCT     50

ref               51 GGCCCTGACCCAGACCTGGGCGGGCTCCCACTCCATGAGGTATTTCTTCA    100
                     ||||||||||||||||||||||||||||||||||||||||||||||||||
seq1              51 GGCCCTGACCCAGACCTGGGCGGGCTCCCACTCCATGAGGTATTTCTTCA    100

ref              101 CATCCGTGTCCCGGCCCGGCCGCGGGGAGCCCCGCTTCATCGCCGTGGGC    150
                     ||||||||||||||||||||||||||||||||||||||||||||||||||
seq1             101 CATCCGTGTCCCGGCCCGGCCGCGGGGAGCCCCGCTTCATCGCCGTGGGC    150

ref              151 TACGTGGACGACACGCAGTTCGTGCGGTTCGACAGCGACGCCGCGAGCCA    200
                     ||||||||||||||||||||||||||||||||||||||||||||||||||
seq1             151 TACGTGGACGACACGCAGTTCGTGCGGTTCGACAGCGACGCCGCGAGCCA    200

ref              201 GAAGATGGAGCCGCGGGCGCCGTGGATAGAGCAGGAGGGGCCGGAGTATT    250
                     ||||||||||||||||||||||||||||||||||||||||||||||||||
seq1             201 GAAGATGGAGCCGCGGGCGCCGTGGATAGAGCAGGAGGGGCCGGAGTATT    250

ref              251 GGGACCAGGAGACACGGAATATGAAGGCCCACTCACAGACTGACCGAGCG    300
                     ||||||||||||||||||||||||||||||||||||||||||||||||||
seq1             251 GGGACCAGGAGACACGGAATATGAAGGCCCACTCACAGACTGACCGAGCG    300

ref              301 AACCTGGGGACCCTGCGCGGCTACTACAACCAGAGCGAGGACGGTTCTCA    350
                     ||||||||||||||||||||||||||||||||||||||||||||||||||
seq1             301 AACCTGGGGACCCTGCGCGGCTACTACAACCAGAGCGAGGACGGTTCTCA    350

ref              351 CACCATCCAGATAATGTATGGCTGCGACGTGGGGCCGGACGGGCGCTTCC    400
                     ||||||||||||||||||||||||||||||||||||||||||||||||||
seq1             351 CACCATCCAGATAATGTATGGCTGCGACGTGGGGCCGGACGGGCGCTTCC    400

ref              401 TCCGCGGGTACCGGCAGGACGCCTACGACGGCAAGGATTACATCGCCCTG    450
                     ||||||||||||||||||||||||||||||||||||||||||||||||||
seq1             401 TCCGCGGGTACCGGCAGGACGCCTACGACGGCAAGGATTACATCGCCCTG    450

ref              451 AACGAGGACCTGCGCTCTTGGACCGCGGCGGACATGGCAGCTCAGATCAC    500
                     ||||||||||||||||||||||||||||||||||||||||||||||||||
seq1             451 AACGAGGACCTGCGCTCTTGGACCGCGGCGGACATGGCAGCTCAGATCAC    500

ref              501 CAAGCGCAAGTGGGAGGCGGTCCATGCGGCGG-AGCAGCGGAGAGTCTAC    549
                     |||||||||||||||||||||||||||||||| |||||||||||||||||
seq1             501 CAAGCGCAAGTGGGAGGCGGTCCATGCGGCGGGAGCAGCGGAGAGTCTAC    550

ref              550 CTGGAGGGCCGGTGCGTGGACGGGCTCCGCAGATACCTGGAGAACGGGAA    599
                     ||||||||||||||||||||||||||||||||||||||||||||||||||
seq1             551 CTGGAGGGCCGGTGCGTGGACGGGCTCCGCAGATACCTGGAGAACGGGAA    600

#---------------------------------------
#---------------------------------------
> writePairwiseAlignments(res2)
########################################
# Program: Biostrings (version 2.38.4), a Bioconductor package
# Rundate: Thu Dec 14 11:35:28 2017
########################################
#=======================================
#
# Aligned_sequences: 2
# 1: ref
# 2: seq2
# Matrix: NA
# Gap_penalty: 10.5
# Extend_penalty: 0.5
#
# Length: 1121
# Identity:     490/1121 (43.7%)
# Similarity:    NA/1121 (NA%)
# Gaps:         577/1121 (51.5%)
# Score: 1882.5
#
#
#=======================================

ref                1 ATGGCCGTCATGGCGCCCCGAACCCTCCTCCTGCTACTCTCGGGGGCCCT     50
                                                                       
seq2               1 --------------------------------------------------      0

ref               51 GGCCCTGACCCAGACCTGGGCGGGCTCCCACTCCATGAGGTATTTCTTCA    100
                                            |||||||||||||||||||||||  ||
seq2               1 -----------------------GCTCCCACTCCATGAGGTATTTCCACA     27

ref              101 CATCCGTGTCCCGGCCCGGCCGCGGGGAGCCCCGCTTCATCGCCGTGGGC    150
                     |  || ||||||||||||||||||||||||||||||||||| ||||||||
seq2              28 CCGCCATGTCCCGGCCCGGCCGCGGGGAGCCCCGCTTCATCACCGTGGGC     77

ref              151 TACGTGGACGACACGCAGTTCGTGCGGTTCGACAGCGACGCCGCGAGCCA    200
                     |||||||||||||||| ||||||| ||||||||||||||||| |||| |
seq2              78 TACGTGGACGACACGCTGTTCGTGAGGTTCGACAGCGACGCCACGAGTCC    127

ref              201 GAAGATGGAGCCGCGGGCGCCGTGGATAGAGCAGGAGGGGCCGGAGTATT    250
                     || || ||||||||||||||| ||||||||||||||||||||||||||||
seq2             128 GAGGAAGGAGCCGCGGGCGCCATGGATAGAGCAGGAGGGGCCGGAGTATT    177

ref              251 GGGACCAGGAGACACGGA------------------ATATG---AAGGCC    279
                     |||||| |||||||| ||                  | ||    ||| ||
seq2             178 GGGACCGGGAGACACAGATCTCCAAGACCGAGACACAGATCTCCAAGACC    227

ref              280 CACTCACAGACTGACCGAGCGAACCTGGGGACCCTGCGCGGCTACTACAA    329
                      || |||||||| |||||| || |||| ||| ||||||||||||||||||
seq2             228 AACACACAGACTTACCGAGAGAGCCTGCGGAACCTGCGCGGCTACTACAA    277

ref              330 CCAGAGCGAGGACGGTTCTCACACCATCCAGATAATGTATGGCTGCGACG    379
                     ||||||||||| ||| ||||||||| ||||||  ||||| ||||||||||
seq2             278 CCAGAGCGAGGCCGGGTCTCACACCCTCCAGAGCATGTACGGCTGCGACG    327

ref              380 TGGGGCCGGACGGGCGCTTCCTCCGCGGG--TACCGGCAGGACGCCTACG    427
                     ||||||||||||||||| |||||||||||  || |  ||| |||||||||
seq2             328 TGGGGCCGGACGGGCGCCTCCTCCGCGGGCATAAC--CAGTACGCCTACG    375

ref              428 ACGGCAAGGATTACATCGCCCTGAACGAGGACCTGCGCTCTTGGACCGCG    477
                     |||||||||||||||||||||||||||||||||||||||| |||||||||
seq2             376 ACGGCAAGGATTACATCGCCCTGAACGAGGACCTGCGCTCCTGGACCGCG    425

ref              478 GCGGACATGGCAGCTCAGATCACCAAGCGCAAGTGGGAGGCGGTCCATGC    527
                     |||||||  || |||||||||||| |||||||||||||||||| || ||
seq2             426 GCGGACACCGCGGCTCAGATCACCCAGCGCAAGTGGGAGGCGGCCCGTGT    475

ref              528 GGCGGAGCAGCGGAGAGTCTACCTGGAGGGCCGGTGCGTGGACGGGCTCC    577
                     ||||||||||   |||| |||||||||||||  |||||||||  ||||||
seq2             476 GGCGGAGCAGGACAGAGCCTACCTGGAGGGCACGTGCGTGGAGTGGCTCC    525

ref              578 GCAGATACCTGGAGAACGGGAAGGAGACGCTGCAGCGCACGGACCCCCCC    627
                     ||||||||||||||||||||||||| |||||| ||||| |||        
seq2             526 GCAGATACCTGGAGAACGGGAAGGACACGCTGGAGCGCGCGG--------    567


ref             1078 CTCACAGCTTGTAAAGTGTGA   1098
                                          
seq2             568 ---------------------    567

#---------------------------------------
#---------------------------------------
>

>ref
ATGGCCGTCATGGCGCCCCGAACCCTCCTCCTGCTACTCTCGGGGGCCCTGGCCCTGACC
CAGACCTGGGCGGGCTCCCACTCCATGAGGTATTTCTTCACATCCGTGTCCCGGCCCGGC
CGCGGGGAGCCCCGCTTCATCGCCGTGGGCTACGTGGACGACACGCAGTTCGTGCGGTTC
GACAGCGACGCCGCGAGCCAGAAGATGGAGCCGCGGGCGCCGTGGATAGAGCAGGAGGGG
CCGGAGTATTGGGACCAGGAGACACGGAATATGAAGGCCCACTCACAGACTGACCGAGCG
AACCTGGGGACCCTGCGCGGCTACTACAACCAGAGCGAGGACGGTTCTCACACCATCCAG
ATAATGTATGGCTGCGACGTGGGGCCGGACGGGCGCTTCCTCCGCGGGTACCGGCAGGAC
GCCTACGACGGCAAGGATTACATCGCCCTGAACGAGGACCTGCGCTCTTGGACCGCGGCG
GACATGGCAGCTCAGATCACCAAGCGCAAGTGGGAGGCGGTCCATGCGGCGGAGCAGCGG
AGAGTCTACCTGGAGGGCCGGTGCGTGGACGGGCTCCGCAGATACCTGGAGAACGGGAAG
GAGACGCTGCAGCGCACGGACCCCCCCAAGACACATATGACCCACCACCCCATCTCTGAC
CATGAGGCCACCCTGAGGTGCTGGGCCCTGGGCTTCTACCCTGCGGAGATCACACTGACC
TGGCAGCGGGATGGGGAGGACCAGACCCAGGACACGGAGCTCGTGGAGACCAGGCCTGCA
GGGGATGGAACCTTCCAGAAGTGGGCGGCTGTGGTGGTGCCTTCTGGAGAGGAGCAGAGA
TACACCTGCCATGTGCAGCATGAGGGTCTGCCCAAGCCCCTCACCCTGAGATGGGAGCTG
TCTTCCCAGCCCACCATCCCCATCGTGGGCATCATTGCTGGCCTGGTTCTCCTTGGAGCT
GTGATCACTGGAGCTGTGGTCGCTGCCGTGATGTGGAGGAGGAAGAGCTCAGATAGAAAA
GGAGGGAGTTACACTCAGGCTGCAAGCAGTGACAGTGCCCAGGGCTCTGATGTGTCTCTC
ACAGCTTGTAAAGTGTGA
>seq1
ATGGCCGTCATGGCGCCCCGAACCCTCCTCCTGCTACTCTCGGGGGCCCTGGCCCTGACC
CAGACCTGGGCGGGCTCCCACTCCATGAGGTATTTCTTCACATCCGTGTCCCGGCCCGGC
CGCGGGGAGCCCCGCTTCATCGCCGTGGGCTACGTGGACGACACGCAGTTCGTGCGGTTC
GACAGCGACGCCGCGAGCCAGAAGATGGAGCCGCGGGCGCCGTGGATAGAGCAGGAGGGG
CCGGAGTATTGGGACCAGGAGACACGGAATATGAAGGCCCACTCACAGACTGACCGAGCG
AACCTGGGGACCCTGCGCGGCTACTACAACCAGAGCGAGGACGGTTCTCACACCATCCAG
ATAATGTATGGCTGCGACGTGGGGCCGGACGGGCGCTTCCTCCGCGGGTACCGGCAGGAC
GCCTACGACGGCAAGGATTACATCGCCCTGAACGAGGACCTGCGCTCTTGGACCGCGGCG
GACATGGCAGCTCAGATCACCAAGCGCAAGTGGGAGGCGGTCCATGCGGCGGGAGCAGCG
GAGAGTCTACCTGGAGGGCCGGTGCGTGGACGGGCTCCGCAGATACCTGGAGAACGGGAA
GGAGACGCTGCAGCGCACGGACCCCCCCAAGACACATATGACCCACCACCCCATCTCTGA
CCATGAGGCCACCCTGAGGTGCTGGGCCCTGGGCTTCTACCCTGCGGAGATCACACTGAC
CTGGCAGCGGGATGGGGAGGACCAGACCCAGGACACGGAGCTCGTGGAGACCAGGCCTGC
AGGGGATGGAACCTTCCAGAAGTGGGCGGCTGTGGTGGTGCCTTCTGGAGAGGAGCAGAG
ATACACCTGCCATGTGCAGCATGAGGGTCTGCCCAAGCCCCTCACCCTGAGATGGGAGCT
GTCTTCCCAGCCCACCATCCCCATCGTGGGCATCATTGCTGGCCTGGTTCTCCTTGGAGC
TGTGATCACTGGAGCTGTGGTCGCTGCCGTGATGTGGAGGAGGAAGAGCTCAGATAGAAA
AGGAGGGAGTTACACTCAGGCTGCAAGCAGTGACAGTGCCCAGGGCTCTGATGTGTCTCT
CACAGCTTGTAAAGTGTGA
>seq2
GCTCCCACTCCATGAGGTATTTCCACACCGCCATGTCCCGGCCCGGCCGCGGGGAGCCCC
GCTTCATCACCGTGGGCTACGTGGACGACACGCTGTTCGTGAGGTTCGACAGCGACGCCA
CGAGTCCGAGGAAGGAGCCGCGGGCGCCATGGATAGAGCAGGAGGGGCCGGAGTATTGGG
ACCGGGAGACACAGATCTCCAAGACCGAGACACAGATCTCCAAGACCAACACACAGACTT
ACCGAGAGAGCCTGCGGAACCTGCGCGGCTACTACAACCAGAGCGAGGCCGGGTCTCACA
CCCTCCAGAGCATGTACGGCTGCGACGTGGGGCCGGACGGGCGCCTCCTCCGCGGGCATA
ACCAGTACGCCTACGACGGCAAGGATTACATCGCCCTGAACGAGGACCTGCGCTCCTGGA
CCGCGGCGGACACCGCGGCTCAGATCACCCAGCGCAAGTGGGAGGCGGCCCGTGTGGCGG
AGCAGGACAGAGCCTACCTGGAGGGCACGTGCGTGGAGTGGCTCCGCAGATACCTGGAGA
ACGGGAAGGACACGCTGGAGCGCGCGG
> sessionInfo()
R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.1 LTS

locale:
 [1] LC_CTYPE=en_IN.UTF-8      
 [2] LC_NUMERIC=C              
 [3] LC_TIME=en_IN.UTF-8       
 [4] LC_COLLATE=en_IN.UTF-8    
 [5] LC_MONETARY=en_IN.UTF-8   
 [6] LC_MESSAGES=en_IN.UTF-8   
 [7] LC_PAPER=en_IN.UTF-8      
 [8] LC_NAME=C                 
 [9] LC_ADDRESS=C              
[10] LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_IN.UTF-8
[12] LC_IDENTIFICATION=C       

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices
[6] utils     datasets  methods   base     

other attached packages:
[1] Biostrings_2.38.4   XVector_0.10.0     
[3] IRanges_2.4.8       S4Vectors_0.8.11   
[5] BiocGenerics_0.16.1

loaded via a namespace (and not attached):
[1] zlibbioc_1.16.0      BiocInstaller_1.20.3
[3] tools_3.2.3

 

1) For seq1 Parameters are kept exactly same as in needle of emboss, still the gaps introduced is at different position.

2) For seq2 even after using the insertion() and deletion() , the reported bases overlap with each other

PS : I have deleted last section of the alignment that is printed because of the word limit for this post.

 
biostrings alignment pairwisealignment • 2.0k views
ADD COMMENT •
0
Entering edit mode
@herve-pages-1542
Last seen 1 day ago
Seattle, WA, United States

Hi Dhwani,

1) One possible explanation for the discrepancy with Emboss needle is that pairwiseAlignment() interprets gapOpening as the cost of opening the gap only. This differs from other alignment software where it's the cost of the first dash ("-") in the gap. Have a look at writePairwiseAlignments() man page (?writePairwiseAlignments) for a detailed explanation of the meaning of gapOpening and the relationship between gapOpening/gapExtension and Gap_penalty/Extend_penalty. Also example C in the Examples section is taken from the Emboss website. It illustrates this relationship and shows you how you need to adjust gapOpening in order to obtain the same results as with Emboss needle.

2) The ranges are reported with respect to the original (i.e. unaligned) pattern. As a consequence, they sometimes seem to overlap:

library(Biostrings)
pa <- pairwiseAlignment("AAACCAAA", "AAATTTTCCTTTTAAA")

insertion(pa)
#  IRangesList of length 1
#  [[1]]
#  IRanges object with 0 ranges and 0 metadata columns:
#         start       end     width
#     <integer> <integer> <integer>

deletion(pa)
# IRangesList of length 1
# [[1]]
# IRanges object with 2 ranges and 0 metadata columns:
#          start       end     width
#      <integer> <integer> <integer>
#  [1]         4         7         4
#  [2]         6         9         4

This is normal and expected. I think that what's confusing here is to think of the start/end/width triplet as representing a range. It makes more sense to interpret the start as the position where the insertion or deletion starts in the original pattern, and the width as its length. The end should be ignored since it has no meaning.  

> writePairwiseAlignments(pa)
########################################
# Program: Biostrings (version 2.47.0), a Bioconductor package
# Rundate: Wed Dec 13 08:25:18 2017
########################################
#=======================================
#
# Aligned_sequences: 2
# 1: P1
# 2: S1
# Matrix: NA
# Gap_penalty: 14.0
# Extend_penalty: 4.0
#
# Length: 16
# Identity:       8/16 (50.0%)
# Similarity:    NA/16 (NA%)
# Gaps:           8/16 (50.0%)
# Score: -39.46618
#
#
#=======================================

P1                 1 AAA----CC----AAA      8
                     |||    ||    |||
S1                 1 AAATTTTCCTTTTAAA     16


#---------------------------------------
#---------------------------------------

Hope this helps.

Cheers,

H.
ADD COMMENT •
0
Entering edit mode

I am able to resolve the confusion for the second query with the help of your explanations. However i tried addressing the first problem based on the formula for gap opening and extension. The answer still remains same even after adding the changes. I have used Gap opening penalty as 10.5 (10 opening + .5 extension) and .5 as gap extension penalty. (While using needle of emboss the gap opening and extension penalty are 10 and .5 respectively. )

ADD REPLY •
0
Entering edit mode

Did you consult ?writePairwiseAlignments as I suggested? It explains the relationship between Biostrings' gapOpening/gapExtension and Emboss needle's Gap_penalty/Extend_penalty, and gives the formula to go from the latter to the former:

    Gap_penalty = gapOpening + gapExtension
    Extend_penalty = gapExtension

So if you use a Gap_penalty and Extend_penalty of 10 and .5 with Emboss needle, you should use a gapOpening and gapExtension of 9.5 and 0.5 with Biostrings.

ADD REPLY •
0
Entering edit mode

The man page for writePairwiseAlignments even contain the following example:

## --------------------------------------------------------------
## C. REPRODUCING THE ALIGNMENT SHOWN AT
##    http://emboss.sourceforge.net/docs/themes/alnformats/align.pair
## --------------------------------------------------------------
pattern <- c("TSPASIRPPAGPSSRPAMVSSRRTRPSPPGPRRPTGRPCCSAAPRRPQAT",
             "GGWKTCSGTCTTSTSTRHRGRSGWSARTTTAACLRASRKSMRAACSRSAG",
             "SRPNRFAPTLMSSCITSTTGPPAWAGDRSHE")
subject <- c("TSPASIRPPAGPSSRRPSPPGPRRPTGRPCCSAAPRRPQATGGWKTCSGT",
             "CTTSTSTRHRGRSGWRASRKSMRAACSRSAGSRPNRFAPTLMSSCITSTT",
             "GPPAWAGDRSHE")
pattern <- unlist(AAStringSet(pattern))
subject <- unlist(AAStringSet(subject))
pattern  # original pattern
subject  # original subject
data(BLOSUM62)
pa5 <- pairwiseAlignment(pattern, subject,
                         substitutionMatrix=BLOSUM62,
                         gapOpening=9.5, gapExtension=0.5)
pa5
writePairwiseAlignments(pa5, Matrix="BLOSUM62")

You can't say the information is not there ;-)

ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 686 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6