Удалить символ из середины строки
У меня есть файл SAM с полем RX:, содержащим 12 оснований, разделенных посередине -, т.е. RX:Z:CTGTGC-TCGTAA
Я хочу удалить дефис из этого поля, но я не могу просто удалить все дефисы из всего файла, так как они содержатся в прочитанных именах, например 1713704_EP0004-T
В основном пытались tr,, но это просто удаление всех дефисов из файла.:
tr -d '"-' < sample.fq.unaln.umi.sam > sample.fq.unaln.umi.re.sam
input представляет собой большой SAM-файл из >10 000 000 строк, например:
1902336-103-016_C1D1_1E-T:34 99 chr1 131341 36 146M = 131376 182 GGACAGGGAGTGTTGACCCTGGGCGGCCCCCTGGAGCCACCTGCCCTGAAAGCCCAGGGCCCGCAACCCCACACACTTTGGGGCTGGTGGAACCTGGTAAAAGCTCACCTCCCACCATGGAGGAGGAGCCCTGGGCCCCTCAGGGG NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN MC:Z:147M MD:Z:83T62cD:i:4 cE:f:0 PG:Z:bwa RG:Z:A MI:Z:34 NM:i:1 cM:i:3 MQ:i:36 UQ:i:45 AS:i:141 XS:i:136 RX:Z:CTGTGC-TCGTAA
Желаемый результат (т.е. последнее поле)
1902336-103-016_C1D1_1E-T:34 99 chr1 131341 36 146M = 131376 182 GGACAGGGAGTGTTGACCCTGGGCGGCCCCCTGGAGCCACCTGCCCTGAAAGCCCAGGGCCCGCAACCCCACACACTTTGGGGCTGGTGGAACCTGGTAAAAGCTCACCTCCCACCATGGAGGAGGAGCCCTGGGCCCCTCAGGGG NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN MC:Z:147M MD:Z:83T62cD:i:4 cE:f:0 PG:Z:bwa RG:Z:A MI:Z:34 NM:i:1 cM:i:3 MQ:i:36 UQ:i:45 AS:i:141 XS:i:136 RX:Z:CTGTGCTCGTAA
Как решить эту проблему?
В частности, неясно, что такое «прочитанное имя» и где оно встречается в строке ввода.
Пожалуйста, опубликуйте несколько строк из файла и опубликуйте пример вывода.
@sjsam Готово, надеюсь понятно
Это еще не ясно. Пожалуйста, опубликуйте строки, которые требуют замены НЕТ. Соответствует ли последнее поле шаблону?
@sjsam Последнее поле — это то, что я хочу изменить RX:Z:CTGTGC-TCGTAA. Моя проблема в том, что я не могу просто удалить все «-» из файла, потому что в прочитанном имени есть дефисы 1902336-103-016_C1D1_1E-T:34, то есть первое поле
"заголовок"? Вы можете уточнить?
Включая мои три ответа были даны на ваш вопрос, попробуйте их и дайте нам отзывы, пожалуйста.
Жалею, что не уделил много внимания 10,000,000. Это действительно большой файл, который превышает 5190 МБ, учитывая длину вашей записи. Надеюсь, у вас есть SSD и приличный процессор. Сообщите мне результаты.
Вы можете принять ответ, если он помог
Почему у вас вообще есть файл SAM? Вы должны никогда иметь файлы SAM, только файлы BAM. Если у вас есть файл BAM, правильное решение относительно легко использовать, например. htslib и несколько строк кода на C. Ответ, который вы сейчас приняли, напротив, неверен.
@KonradRudolph У меня был файл sam, потому что мне нужно было внести эти изменения. Я признаю, что у меня нет опыта работы с pysam или любым парсером bam, поэтому я преобразовал свой bam в sam, чтобы помочь решить мою проблему. Эти правки на месте работают, поэтому я бы не назвал их «неправильными», просто не оптимальными. Я закончил тем, что применил регулярное выражение к проблеме, которая является гораздо более безопасным способом сделать это. Я сейчас пытаюсь выучить пысам, так как мои знания С или Perl равны 0.
@lgallagher Я назвал принятый ответ неправильным, потому что он заменит -где угодно, а не только в выбранном вами теге. Это очень неправильно.
@KonradRudolph Я обновил свой вопрос с помощью решения pysam. Приятно получать отзывы о нем
@lgallagher, пожалуйста, опубликуйте это как ответ (самостоятельные ответы приветствуются!)
@KonradRudolph Ах, извините - все еще разбираюсь с stackoverflow
Ответы 4
аук
awk '{sub(/-/,"",$NF)}1' file
это то, что вам нужно.
Объяснение
- Из это видно, что вас интересует только последнее поле.
- NF — это общее количество полей, содержащихся в записи, поэтому $NF — это последнее поле.
sub(/-/,"",$NF)заменяет-в последнем поле пустой строкой, делая изменение постоянным.
GNU-сед
По той же причине это,
sed -Ei 's/^(.*)-/\1/' file
буду работать. У него есть дополнительное преимущество, заключающееся в том, что он может выполнять редактирование на месте.
Объяснение
- Параметр
-Eвключает расширенный движок регулярных выражений. (.*)— это жадный поиск, который будет соответствовать любому символу (.) любое количество раз (*). Для того, чтобы быть жадным, он будет соответствовать чему угодно до последнего дефиса.()заставляетsedпомнить, что было сопоставлено.- В части замены мы помещаем только совпадающую часть
\1(1, потому что у нас есть только одна пара скобок, обратите внимание, что вы можете иметь столько, сколько хотите) без дефиса, таким образом эффективно удаляя ее из последнего поля, где она должна встречаться .
Note : The GNU awk support -i inplace, but I'm not sure from which version on.
1 это опечатка?
@PaulHodges: Нет, 1 просто печатает запись. Быть идиоматичным :-)
Интересно. Почему печатает запись? Все еще учусь awk. Я бы просто явно добавил print. :)
@Paul, это сокращение от 1{print}, когда условие оценивается как true/1, его действие выполняется. Если условие не имеет действия, по умолчанию используется {print}.
@PaulHodges: 1 — самая простая команда, которую понимает awk. Он просто печатает запись. Попробуйте awk '1' file :)
Я просто удивился, увидев его в конце блока кода, но тогда я думаю, что технически он находится в начале подразумеваемого {print}, а?
Вы не хотите делать это на всех строках 10-мегабайтного файла.
@Иниан. В этом случае он может попробовать мое sed решение. :-) В любом случае я не утверждаю, что это будет очень быстро. Но это может совпасть с каким-нибудь талармейным python решением. Разве это не так?
Ваше решение sed - это именно то, что мне нужно. Из-за размера моих файлов довольно сложно найти «быстрое» решение. Я могу запустить это на узле HPC, так что это должно помочь. Спасибо тебе за помощь
wrt [sed] has an added advantage that it can perform an inplace edit. — вы используете GNU sed, чтобы получить -i, и вы можете использовать GNU awk, чтобы получить эквивалент -i inplace. Эта конкретная команда awk изменит пробелы в остальной части вывода, что может быть нежелательно.
@EdMorton: Спасибо, обновляю ответ. У меня есть gawk 4.02, в котором нет опции -i в руководстве. Не уверен, когда он был введен.
-i в gawk не означает то же самое, что и в sed. В sed -i означает inplace, но в gawk -i означает include, а inplace — это подключаемая библиотека, поэтому эквивалентом sed -i является awk -i inplace. -i появился в gawk 4.1, я не уверен, появилась ли библиотека inplace в том же выпуске или позже. Сейчас у нас gawk 5.0, так что в любом случае это было давно!
@EdMorton Я немного удивлен, что мой CentOS 7 все еще работает на более старой версии awk. У него также нет опции -i.
Это очень плохо. Если вам интересно посмотреть историю того, когда/почему/как появился -i inplace, все началось в списке рассылки давным-давно на сервере далеко-далеко... lists.gnu.org/archive/html/bug-gawk/2012-12/msg00048.html
Осторожнее с этим ответом, он неправильный: он заменяет дефисы все в последнем столбце. Но в реальных файлах SAM этот столбец содержит несколько полей структурированных данных, и замена должна происходить только в одном из этих полей, а не в любом другом.
@KonradRudolph Из это я понимаю, что оператор не предоставил нам надлежащего образца для работы. awk заменит только первый дефис. Но, как я смиренно признаю, это не правильные инструменты для работы.
@sjsam Да, ОП оставил неявными многие предположения, о которых вы знаете, только если у вас есть опыт в биоинформатике.
Этот шаблон есть во многих записях, которые вы хотите отредактировать, и всегда в конце строки? Если так -
sed -E 's/^(.*)(\s..:.:......)-(......\s*)$/\1\2\3/' < sample.fq.unaln.umi.sam > sample.fq.unaln.umi.re.sam
Он не всегда будет в конце строки, и количество букв до/после тире не будет фиксированным.
Лучшее решение — работать с файлами BAM, а не с файлами SAM, и использовать подходящую библиотеку парсера/записи BAM, такую как htslib.
Если этого не хватает, вы можете собрать что-нибудь вместе, выполнив поиск регулярного выражения ^RX:Z: в необязательных тегах (столбцы 12 и выше).
Работа со столбцами, хотя и возможна, с sed затруднена. Вместо этого вот как это сделать в awk:
awk -F '[[:space:]]*' '{
for (i = 12; i <= NF; i++) {
if ($i ~ /^RX:Z:/) gsub("-", "", $i)
}
}
1' file.sam
А вот примерно эквивалентное решение в виде «однострочника» на Perl:
perl -ape '
for (@F[11..(scalar @F)]) {
s/-//g if /^RX:Z:/;
}
$_ = join("\t", @F);
' file.sam
Чтобы выполнить замену в исходном файле, вы можете передать опцию -i.bak в perl (это создаст резервную копию file.sam.bak; если вы не хотите резервную копию, опустите расширение).
Я решил эту проблему, используя pysam, который быстрее, безопаснее и требует меньше места на диске, так как файл sam не требуется. Это не идеально, я все еще изучаю python и использую pysam полдня.
import pysam
import sys
from re import sub
# Provide a bam file
if len(sys.argv) == 2:
assert sys.argv[1].endswith('.bam')
# Makes output filehandle
inbamfn = sys.argv[1]
outbamfn = sub('.bam$', '.fixRX.bam', inbamfn)
inbam = pysam.Samfile(inbamfn, 'rb')
outbam = pysam.Samfile(outbamfn, 'wb', template=inbam)
# Counters for reads processed and written
n = 0
w = 0
# .get_tag() retrieves RX tag from each read
for read in inbam.fetch(until_eof=True):
n += 1
umi = read.get_tag('RX')
assert umi is not None
umifix = umi[:6] + umi[7:]
read.set_tag('RX', umifix, value_type='Z')
if '-' in umifix:
print('Hyphen found in UMI:', umifix, read)
break
else:
w += 1
outbam.write(read)
inbam.close()
outbam.close()
print ('Processed', n, 'reads:\n',
w, 'UMIs written.\n',
str(int((w / n) * 100)) + '% of UMIs fixed')
Инструмент для работы - это всегда прелесть!
Размещение некоторых образцов строк было бы полезно.