Я занимаюсь геномикой. У меня файл с форматом FASTA читает. Это гены. Каждый ген называется считыванием или контигом. Каждый контиг начинается с заголовка и сопровождается алфавитами или нуслеотидами, например ACTG, определенной длины. Я хочу определить самый длинный контиг и самый короткий контиг, чтение или ген в этом файле. Подскажите пожалуйста сценарий ubuntu, чтобы найти такие контиги. Каждый контиг или чтение в этом формате FASTA выглядит следующим образом:
>Locus_1000_Transcript_1/1_Confidence_0.000_Length_648 FTBs=645 (Header)
CcGccttggtaacctCgccAGcatATtgagcTttGGatccGGaTggtcgtaGaAtgGCaaG
GcaGgagAgAgtgtctaatgtggCgccGctctgtAccCgGgGGgTAACaAtgAATTtGCga
CgaCGtggTAtGcCcttCGttgAaacccTtaTtagttGgAGCcGctAtgtggcgGTccaat
TaTcaagtAttTcCCACaTcttgAagCgcttcTgGATgTacgCatactatgggTtgacgtt
AGtGtAgCcgAgattTCacaGtAgctcCGAACGgtgGTagCAgacGcccGttCacAAaAaC
Заголовок имеет определенный формат, который показывает локусы генов и количество генов, и между каждым контигом или прочтением будет пробел. Каждое мое чтение или контиг в файле будет начинаться с заголовка того же типа, что и упомянутый выше, но значения могут отличаться. Каждый контиг или чтение начинается со знака>. Могут быть контиги одинаковой длины. - Наука 3 минуты назад
Принятие Length
значения в заголовках FASTA правильны, я извлек бы их оттуда:
sed -nre 's/^>.*_Length_([0-9]+) .*/\1/p' \
затем отсортируйте их численно
| sort -n \
затем произведите первую и последнюю строку
| sed -ne '1p;$p'
В одном операторе:
sed -nre 's/^>.*Length_([0-9]+) .*/\1/p' | sort -n | sed -ne '1p;$p'
Если длинам, объявленным в заголовках, нельзя доверять, то считать длину последовательностей FASTA, я сначала преобразовал бы их в unfasta, то распечатал бы длину строки каждой второй строки в то же sort | sed
отфильтруйте как выше:
uf | awk 'NR%2==0 {print length}' | sort -n | sed -n '1p;$p'
Где uf
простой сценарий удара, найденный здесь.
Примечание: обе остроты являются фильтрами, который является, они читают свой вход из стандартного входа и пишут в стандартный вывод. Использовать cat
подавать их файлы (или wget -O -
подавать их от Интернета).
Этот сценарий Python создает словарь записей и использует линейный поиск для нахождения, какой является самым коротким, какой является самым длинным в файле. Обратите внимание на то, что это игнорирует регистры, где существует два контига того же значения (хотя это может также быть реализовано.
Код:
#!/usr/bin/env python3
import sys
def main():
records = {}
current_length = 0
current_contig = ''
with open(sys.argv[1]) as f:
for index,line in enumerate(f,1):
if line == '\n': continue
if line.startswith('>'):
if current_contig != line:
records[current_contig] = current_length
current_contig = line.strip()
current_length = 0
else:
current_length = current_length + len(line.strip())
records[current_contig] = current_length
records.pop('')
shortest_contig = None
longest_contig = None
longest_val = 0
shortest_val = float("inf")
for contig,length in records.items():
if length < shortest_val:
shortest_val = length
shortest_contig = contig
if length > longest_val:
longest_val = length
longest_contig = contig
print('Longest: ' + longest_contig)
print('Shortest: ' + shortest_contig)
if __name__ == '__main__': main()
Тестовый прогон:
$ cat input.txt
> Entry 1
CcGccttggtaacctCgccAGcatATtgagcTttGGatccGGaTggtcgtaGaAtgGCaaG
GcaGgagAgAgtgtctaatgtggCgccGctctgtAccCgGgGGgTAACaAtgAATTtGCga
CgaCGtggTAtGcCcttCGttgAaacccTtaTtagttGgAGCcGctAtgtggcgGTccaat
TaTcaagtAttTcCCACaTcttgAagCgcttcTgGATgTacgCatactatgggTtgacgtt
AGtGtAgCcgAgattTCacaGtAgctcCGAACGgtgGTagCAgacGcccGttCacAAaAaC
> Entry 2
CcGccttggtaacctCgccAGcatATtgagcTttGGatccGGaTggtcgtaGaAtgGCaaG
GcaGgagAgAgtgtctaatgtggCgccGctctgtAccCgGgGGgTAACaAtgAATTtGCga
CgaCGtggTAtGcCcttCGttgAaacccTtaTtagttGgAGCcGctAtgtggcgGTccaat
TaTcaagtAttTcCCACaTcttgAagCgcttcTgGATgTacgCatactatgggTtgacgtt
AGtGtAgCcgAgattTCacaGtAgctcCGAACGgtgGTagCA
$ python3 contigs.py input.txt
Longest: > Entry 1
Shortest: > Entry 2