如何构建fasta文件及其index索引文件

需要准备的文件是序列信息bed格式文件或者GRanges格式文件

比如有一个bed文件：

$ cat test.bed
chr1    1000010 1000020
chr2    1000021 1000030
chr3    1000031 1000040

我习惯保存为GR格式：

> test = read.table("test.bed")
> gr = toGR(test)
> gr
GRanges object with 3 ranges and 0 metadata columns:
      seqnames          ranges strand
         <Rle>       <IRanges>  <Rle>
  [1]     chr1 1000010-1000020      *
  [2]     chr2 1000021-1000030      *
  [3]     chr3 1000031-1000040      *
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

toGR是我自己编写的函数，这里就不放出来了。此外还编写了另外两个函数，分别为：

getSeq，用于得到碱基序列

getSeq <- function(POSs){
    
    suppressPackageStartupMessages(library(BSgenome.Hsapiens.UCSC.hg19))
    seqs = c()
    for(i in 1:length(POSs)){
        POS = POSs[i]
        chr = sapply(strsplit(POS, "[-:]"), function(z) z[1])
        strat = as.numeric(sapply(strsplit(POS, "[-:]"), function(z) z[2]))
        end = as.numeric(sapply(strsplit(POS, "[-:]"), function(z) z[3]))
        seqs = c(seqs, as.character(BSgenome.Hsapiens.UCSC.hg19[[chr]][strat:end]))
    }
    return(unlist(seqs))
}

printSeq，用于打印碱基序列

printSeq <- function(Seq, binwidth = 50){
   N_line = nchar(Seq) / binwidth 
   if(N_line > as.integer(N_line)){
       N =  as.integer(N_line) + 1
   } else{
       N =  as.integer(N_line)
   }

   bases = strsplit(Seq, "")[[1]]
   for(i in 0:(N-1)){
       xseq = paste0(bases[(1 + binwidth*i):(binwidth + i*binwidth)], collapse = "") 
       if(i != (N-1)){
           xseq = paste0(bases[(1 + binwidth*i):(binwidth + i*binwidth)], collapse = "") 
       } else{
           xseq = paste0(bases[(1 + binwidth*i):length(bases)], collapse = "") 
       }
       cat(xseq, "\n")
   }
}

默认在fasta中每行显示50个碱基

然后使用前面的GR数据生成fasta文件：

> seqs = getSeq(as.character(gr))
> seqs
[1] "CTCACCCAGGA" "AAATTGAAGA"  "AGAAGGAAAC"
> ID = as.character(seqnames(gr))
> ID
[1] "chr1" "chr2" "chr3"

sink("test.fa")
for(i in 1:length(seqs)){
    cat(paste0(">", ID[i]), "\n")
    printSeq(seqs[i])
}
system("sed 's/ //g' test.fa > tmp.fa")
system("rm test.fa")
system("mv tmp.fa test.fa")
sink()

这样就生成了fasta文件了：

$ cat test.fa 
>chr1
CTCACCCAGGA
>chr2
AAATTGAAGA
>chr3
AGAAGGAAAC

最后还需要创建索引：

system("samtools faidx test.fa")

这样就完成了~

好家伙，samtools有提取序列的功能，还能指定每行的碱基数，位点格式为：chr:start-end，例如

$ samtools faidx hg19.fa chr1:2066354-2066580
>chr1:2066354-2066580
AGGCAGGGGACAGACGGACCCGGCCTGCGTTGGCCTGGGGTGACTTCACGGCTCCACTGT
CAGCAAGCGGCCGTCCCGTGGTGGATCCTGTCCGCCCTGCGAGGACACCTGGCTCCATCC
ACACCTGGGCCTCTGTCTCCAGCCGCCGAGGCCGTGACACCATGAGGATCATGTGAGGAG
GGGCAGAGAGAGGCCTCCGGGAGGCCGTCATTCCAGCCCTGCCTTCC

还能将想要提取的区域写入文件中，并传递给-r参数。每个区域单独一行

区域的格式刚好是Granges对象非常容易生成的格式，所以可以写个GR转fasta的小函数：

grToFasta <- function(gr, file){

    hg19.fa = "hg19.fa"
    pos = as.character(gr)
    write.table(pos, "tmp.pos", sep = "\n", quote = FALSE, col.names = FALSE, row.names = FALSE)
    system(paste("samtools faidx -n 50", hg19.fa, "-r tmp.pos >", file))
    system(paste("rm tmp.pos & samtools faidx", file))
}

我直接好家伙，还是调包香~

最后编辑于：2021.01.03 16:46:52