1.生成一組示例序列
先生成一條序列,並且把生成序列的代碼寫成個函數,函數的輸入數據是序列的長度:
x= c("A","T","C","G")
ms paste(sample(x,n,replace = TRUE),collapse = "")
}
ms(25)
#> [1] "AAAGGTAGGTGTTGTAGAACTATGT"
然後生成一組示例數據,例子是10個序列,每個長度25:
知識點一:用sample和paste0生成隨機字符串
paste0和sample都是基礎函數,但是用處很大,sample表示抽樣,加上replace = TRUE就是有放回的抽樣,paste0則表示將結果無縫連接起來,這不就成了個序列嘛~
知識點二:
用數據結構(例如向量/列表)儲存for循環多次產生的結果,這些結果成為向量/列表裡的元素。
y = c()
for(i in 1:10){
y[[i]] = ms(25)
}
y
#> [1] "GACCGCCTTATAGCAGTTAGTGTAG" "GGAAGCGTCATTCGAACTATGCCGG"
#> [3] "AACAGTTCCGTCCACTTGGCCATTG" "TACCGCCTTTGTTAGCCGTGTGGCT"
#> [5] "ATCTCCAGACTTTATAAATCTATCA" "GTGTGGTATTATCGTGATCTCCACG"
#> [7] "GAGCACAGACGCGCTAGAGAATTGA" "GTTTTACACGCATCAGTGGGTAAAG"
#> [9] "TGGCAGTGTCTCTCTAGCGTGATCT" "AGGCCATGGGTTGGGAGTTACTTAT"
如果你有真實的序列,直接用readLines()讀進來,as.character()變成字符串就好。
2.把序列成fasta格式
最簡單的fasta格式是第一行大於號加序列名稱,第二行是序列。
把字符串序列變成fasta格式,可以這樣
方法一
x1 = paste0("seq",1:10)
res ",x1,"\\n",y)
res
#> [1] ">seq1\\nGACCGCCTTATAGCAGTTAGTGTAG"
#> [2] ">seq2\\nGGAAGCGTCATTCGAACTATGCCGG"
#> [3] ">seq3\\nAACAGTTCCGTCCACTTGGCCATTG"
#> [4] ">seq4\\nTACCGCCTTTGTTAGCCGTGTGGCT"
#> [5] ">seq5\\nATCTCCAGACTTTATAAATCTATCA"
#> [6] ">seq6\\nGTGTGGTATTATCGTGATCTCCACG"
#> [7] ">seq7\\nGAGCACAGACGCGCTAGAGAATTGA"
#> [8] ">seq8\\nGTTTTACACGCATCAGTGGGTAAAG"
#> [9] ">seq9\\nTGGCAGTGTCTCTCTAGCGTGATCT"
#> [10] ">seq10\\nAGGCCATGGGTTGGGAGTTACTTAT"
writeLines(res)
#> >seq1
#> GACCGCCTTATAGCAGTTAGTGTAG
#> >seq2
#> GGAAGCGTCATTCGAACTATGCCGG
#> >seq3
#> AACAGTTCCGTCCACTTGGCCATTG
#> >seq4
#> TACCGCCTTTGTTAGCCGTGTGGCT
#> >seq5
#> ATCTCCAGACTTTATAAATCTATCA
#> >seq6
#> GTGTGGTATTATCGTGATCTCCACG
#> >seq7
#> GAGCACAGACGCGCTAGAGAATTGA
#> >seq8
#> GTTTTACACGCATCAGTGGGTAAAG
#> >seq9
#> TGGCAGTGTCTCTCTAGCGTGATCT
#> >seq10
#> AGGCCATGGGTTGGGAGTTACTTAT
知識點三:
writeLines可以將字符串中的特殊符號"\\n","\\t"等顯示出來成為它本來的樣子。
方法二
(其實寫這一篇只是為了舉個例子說說apply(),結果跑偏了。)
df res2 = apply(df,1,function(x){
paste0(">",x[1],"\\n",x[2])
})
res2
#> [1] ">seq1\\nGACCGCCTTATAGCAGTTAGTGTAG"
#> [2] ">seq2\\nGGAAGCGTCATTCGAACTATGCCGG"
#> [3] ">seq3\\nAACAGTTCCGTCCACTTGGCCATTG"
#> [4] ">seq4\\nTACCGCCTTTGTTAGCCGTGTGGCT"
#> [5] ">seq5\\nATCTCCAGACTTTATAAATCTATCA"
#> [6] ">seq6\\nGTGTGGTATTATCGTGATCTCCACG"
#> [7] ">seq7\\nGAGCACAGACGCGCTAGAGAATTGA"
#> [8] ">seq8\\nGTTTTACACGCATCAGTGGGTAAAG"
#> [9] ">seq9\\nTGGCAGTGTCTCTCTAGCGTGATCT"
#> [10] ">seq10\\nAGGCCATGGGTTGGGAGTTACTTAT"
writeLines(res2)
#> >seq1
#> GACCGCCTTATAGCAGTTAGTGTAG
#> >seq2
#> GGAAGCGTCATTCGAACTATGCCGG
#> >seq3
#> AACAGTTCCGTCCACTTGGCCATTG
#> >seq4
#> TACCGCCTTTGTTAGCCGTGTGGCT
#> >seq5
#> ATCTCCAGACTTTATAAATCTATCA
#> >seq6
#> GTGTGGTATTATCGTGATCTCCACG
#> >seq7
#> GAGCACAGACGCGCTAGAGAATTGA
#> >seq8
#> GTTTTACACGCATCAGTGGGTAAAG
#> >seq9
#> TGGCAGTGTCTCTCTAGCGTGATCT
#> >seq10
#> AGGCCATGGGTTGGGAGTTACTTAT
identical(res,res2)
#> [1] TRUE
知識點四:
apply(X, MARGIN, FUN, …)其中X是數據框/矩陣名;MARGIN為1表示取行,為2表示取列,FUN是函數。強大的apply配上自定義函數,可以做很多事情。
判斷兩個變量是否完全相等,用identical(),不用 ==
3.後續
可以合併
如果你需要合成一整個字符串,那也簡單。
resl=paste(res,collapse = "\\n")
resl
# [1] ">seq1\\nCCTGTCTAGTCCGAGATATGGTAAG\\n>seq2\\nATGTGAAGGAACACGCATGAAAAAG\\n>seq3\\nAGGGCTATCCGGGGAAGAAATTAGC\\n>seq4\\nACGCTATATCGCAGCATGTCTCTAG\\n>seq5\\nATCAGCCACGTCACAACGGTGGCAT\\n>seq6\\nACTACTACTGCTAAAGTAGGCTGCT\\n>seq7\\nGGTTCCATTGTATGACATAACAAAC\\n>seq8\\nTGACAGAACGAAGTCATGATACCGA\\n>seq9\\nATTTGCCTGTGTCCTCGAGAGTAGT\\n>seq10\\nCTCTCAGCGACGGGCGGTGAAAGCT"
writeLines(resl)
#> >seq1
#> GACCGCCTTATAGCAGTTAGTGTAG
#> >seq2
#> GGAAGCGTCATTCGAACTATGCCGG
#> >seq3
#> AACAGTTCCGTCCACTTGGCCATTG
#> >seq4
#> TACCGCCTTTGTTAGCCGTGTGGCT
#> >seq5
#> ATCTCCAGACTTTATAAATCTATCA
#> >seq6
#> GTGTGGTATTATCGTGATCTCCACG
#> >seq7
#> GAGCACAGACGCGCTAGAGAATTGA
#> >seq8
#> GTTTTACACGCATCAGTGGGTAAAG
#> >seq9
#> TGGCAGTGTCTCTCTAGCGTGATCT
#> >seq10
#> AGGCCATGGGTTGGGAGTTACTTAT
還可以導出
write.table(resl,
file = "test.fasta",
row.names = F,
quote = F)
知識點五:
字符串也可以導出的。想要讓導出的文件沒引號,就加quote = F,導出文件後綴可以改,改成fasta沒關係的,因為是純文本,後綴不改變本質,只是決定了這文件在window電腦上的默認打開方式而已。
閱讀更多 生信星球 的文章