Разделение данных на основе численных сравнений в строках FASTA

Question

Я пытаюсь отделить свои данные, которые я импортировал из файла FASTA, используя read.fasta, в список всех данных, где жирный шрифт внизу как минимум в 5 раз больше, чем число справа от него. .

"> NP_000005.2 \ t1474 \ tGI: 66932947 \ tNM_000014.4: 114..4538 \ tGeneID: 2 \ tA2M \ talpha-2-предшественник макроглобулина [Homo sapiens]. \ TSAP2: rs200486795;S1411R; MAF = 0,0119 , 0,025,0.0161 "attr (," class ")

Любые предложения приветствуются.

Answer 1

Я не очень знаком с seqinr, и есть что-то странное, когда он читает заголовки fastta. Итак, я показываю пример с пакетом Biostrings ниже:

library(Biostrings)
fa = readAAStringSet("muscle.var.fasta")
head(names(fa),1)
[1] "NP_000005.2\t1474\tGI:66932947\tNM_000014.4:114..4538\tGeneID:2\tA2M\talpha-2-macroglobulin precursor [Homo sapiens].\tSAP1:rs376343602;T1455M;MAF=0.0243,0.0,0.0168" 

Это пример набора данных:

Data <- c("NP_000005.2\t1474\tGI:66932947\tNM_000014.4:114..4538\tGeneID:2\tA2M\talpha-2-macroglobulin precursor [Homo sapiens].\tSAP1:rs376343602;T1455M;MAF=0.0243,0.0,0.0168", 
"NP_000005.2\t1474\tGI:66932947\tNM_000014.4:114..4538\tGeneID:2\tA2M\talpha-2-macroglobulin precursor [Homo sapiens].\tSAP2:rs200486795;S1411R;MAF=0.0119,0.025,0.0161", 
"NP_000005.2\t1474\tGI:66932947\tNM_000014.4:114..4538\tGeneID:2\tA2M\talpha-2-macroglobulin precursor [Homo sapiens].\tSAP3:rs181129451;R1407W;MAF=0.0,0.7992,0.2571", 
"NP_000005.2\t1474\tGI:66932947\tNM_000014.4:114..4538\tGeneID:2\tA2M\talpha-2-macroglobulin precursor [Homo sapiens].\tSAP4:rs374966730;T1359I;MAF=0.0118,0.0,0.0080", 
"NP_000005.2\t1474\tGI:66932947\tNM_000014.4:114..4538\tGeneID:2\tA2M\talpha-2-macroglobulin precursor [Homo sapiens].\tSAP5:rs372632979;V1345M;MAF=0.0118,0.0,0.0080"
)

Сначала мы преобразуем заголовок в табличную форму,

test = read.delim(text=Data,header=F,stringsAsFactors=F)

и мы хотим, чтобы восьмой столбец содержал MAF, и мы используем strsplit, чтобы разделить эту колонку по символам;=,:

strsplit(test[1,8],"[;=,]")
[[1]]
[1] "SAP1:rs376343602" "T1455M"           "MAF"              "0.0243"          
[5] "0.0"              "0.0168" 

Нам понадобится значение от 4-го до 6-го после разделения, а ниже приведен фрагмент кода, построенный на этом strsplit

AF = t(sapply(strsplit(test[,8],"[;=,]"),function(i)as.numeric(i[4:6])))
> head(AF)
       [,1]   [,2]   [,3]
[1,] 0.0243 0.0000 0.0168
[2,] 0.0119 0.0250 0.0161
[3,] 0.0000 0.7992 0.2571
[4,] 0.0118 0.0000 0.0080
[5,] 0.0118 0.0000 0.0080

Это дает 1-е, 2-еи 3-е значение, после "MAF =", для имен от 1 до 5

, что нам нужно в этом примере, это строки 1,4 и 5, мы можем получить это, используя

wh = which(AF[,1]>5*AF[,2])
#gives the header which fulfills this

test[wh,]
Data[wh]

Давайтепримените это ко всему набору данных

test = read.delim(text=names(fa),header=F,stringsAsFactors=F)
AF = t(sapply(strsplit(test[,8],"[;=,]"),function(i)as.numeric(i[4:6])))
wh = which(AF[,1]>5*AF[,2])
# check hits
length(wh)
[1] 30209
# tail hits
tail(names(fa[wh]),2)
[1] "NP_998839.1\t284\tGI:47519616\tNM_213674.1:240..1094\tGeneID:7169\tTPM2\ttropomyosin beta chain isoform Tpm2.1sm/cy [Homo sapiens].\tSAP4:rs146271535;C36S;MAF=0.0116,0.0,0.0077"
[2] "NP_998890.1\t89\tGI:47524167\tNM_213725.1:130..399\tGeneID:6176\tRPLP1\t60S acidic ribosomal protein P1 isoform 2 [Homo sapiens].\tSAP2:rs149295760;N28S;MAF=0.0116,0.0,0.0077"

Вы можете записать файл фаста или имена:

writeLines(names(fa)[wh],"hits.names.txt")
writeXStringSet(fa[wh],"hits.fa")