readAxt() 将小写转换为大写

Posted 2023-03-23

【中文标题】readAxt() 将小写转换为大写

【英文标题】：readAxt() converts lowercase to uppercase

【发布时间】：2022-01-24 06:47:43

【问题描述】：

有没有什么方法可以读取 AXT 文件而不自动将序列转换为大写？

代码链接为：https://github.com/alexander-nash/kurtosis_conservation/blob/master/get_identical_seq_locations.R

getLengthsOfIdenticalSeqs() 此函数比较两个序列并确定匹配序列的长度。 例如： ATCGCGAT TTCGAAAT 输出： 长度为 3 的 TCG 长度为 3 的 AT

但问题在于 readAxt() 函数会自动将小写转换为大写，然后比较错误的序列。

if(species2 != "lepOcu1")
  axts<-lapply(species2, function(x) 
    lel<-dir(paste0("Human-mouse/Human Mouse 2009/axtNet/"), pattern=paste0("chrX", ".*.axt"), full.names=T)
    lel<-lel[!grepl("Exon", lel)]
    lel<-lel[!grepl("broken", lel)]
    tfn<-paste0("Human-mouse/Human Mouse 2009/bigZips/hg19/", species1, ".2bit")
    if(!file.exists(tfn)) tfn<-paste0("Human-mouse/Human Mouse 2009/bigZips/hg19", species1, "/bigZips/", species1, ".2bit")
    qfn<-paste0("Human-mouse/Human Mouse 2009/bigZips/mm10/", species2, ".2bit")
    if(!file.exists(qfn)) qfn<-paste0("Human-mouse/Human Mouse 2009/bigZips/mm10", species2, "/bigZips/", species2, ".2bit")
    out<-readAxt(lel, tAssemblyFn=NULL, qAssemblyFn=NULL)
  )


names(axts)<-species2
print((axts))

此代码输出以下序列： 具有 80740 个对齐对的 Axt： 1 chrX 70345 70614 chr8 35873813 35874094 - 6175 GGTACTGAGGTCCCCTGGGTACTGAGATCTCCTCGGTACTGAAGTCTCCTCGGTGCTGAGGTCGCCTCGGTGCTG...GGTACTGAGGTCGCCTAGGTACTGAGACCTTCTAGGTCCTGAGGT--------CTAGGTACTGAGG-CCTTCTCC GATGCTGAGGTTCCCAGGATGCTGAGGTTCCCAGGATGCTGAGGTTCCCAGGATGCTGAGGTTCCCAGGATGCTG...GATGCTGAGGTTCCCAGGATGCTGAGGTTCCCAGGATGCTGAGGTTCCTCTCCCAGGATGCTGAGGTTCCTCTCC

但原来的序列是（存在小写）： 0 chrX 70345 70614 chr8 35873813 35874094 - 6175

ggtacTGAGGTCCCCTGGGTACTGAGATCTCCTCGGTACTGAAGTCTCCTCGGTGCTGAGGTCGCCTCGGTGCTGAGACCTCCTAGGTATTGAGGTCGCCTCGGTACTGAGGTTGCCTC----------------------------GGTGCTGAGGT-----CGCCACGGTGCTGAGACCTCCTAGATACTGAGG----TCTCCTAGGCACGGAGATCTCCTATGTACAGAGACCTCGTCGGTACTGAGGTCAGCCTAGGTACTGAGACCTTCTAG-- --CTAGGTACTGAGG-CCTTCTCC

GATGCTGAGGTTCCCAGGATGCTGAGGTTCCCAGGATGCTGAGGTTCCCAGGATGCTGAGGTTCCCAGGATGCTGAGGTTCCCAGGATGCTGAGGTTCCCAGGATGCTGAGGTT-CCTCTCCCGGGATGCTGAGGTTCCTCTCCCGGGATGCTGAGGTTCCTCTCCCAGGATGCTGAGGTTCCCAGGATGCTGAGGTTCCTCTCCCAG --------------------------------- GATGCTGAGGTTCCCAGGATGCTGAGGTTCCCAGGATGCTGAGGTTCCTCTCCCAGGATGCTGAGGTTCCTCTCC P>

【问题讨论】：

请给出一个带有输出和预期输出的简短代码示例。

@AndreWildberg 完成。

【参考方案1】：

据我所知，DNAStringSet() 的 getLengthsOfIdenticalSeqs() 是罪魁祸首。而且它不是预期的错误。

详情： DNAString 类是直接的 XString 子类（没有 额外的插槽）。因此描述的所有功能和方法 在 XString 手册页中也可以使用 DNAString 对象 （继承）。 与 BString 容器不同，它允许存储任何单个 字符串（基于单字节字符集） DNAString 容器只能存储基于 DNA 字母表的字符串（参见 以下）。此外，存储在 DNAString 对象中的字母是 以优化快速搜索算法的方式编码。

您可以将 get_identical_seq_locations.R 中的函数更改为使用readBStringSet() 或尝试在打印输出后将字符串转换回来。

对于后者，从输出中获取序列位置信息很重要。 您可以先保存小写碱基的位置，然后在打印输出后将它们放回原处，例如

DNA <- c("ACGTggTTAa")
lower <- sapply( strsplit( DNA, ""), function(x) grepl("[[:lower:]]",x) )
lower
       [,1]
 [1,] FALSE
 [2,] FALSE
 [3,] FALSE
 [4,] FALSE
 [5,]  TRUE
 [6,]  TRUE
 [7,] FALSE
 [8,] FALSE
 [9,] FALSE
[10,]  TRUE

DNA_out
[1] "ACGTGGTTAA"

DNA_out_split <- unlist(strsplit( DNA_out, "" ))

DNA_out_split[lower] <- tolower( DNA_out_split[lower] )

DNA_out <- paste(DNA_out_split, collapse="")
DNA_out
[1] "ACGTggTTAa"

我只是不知道打印出来后是否有办法知道基本位置。