Введение

Во время выполнения анализа с помощью программного обеспечения Genome Analyzer Pipeline создается несколько промежуточных файлов. В этом примере вы узнаете, как читать и манипулировать информацией, содержащейся в файлах последовательности (_sequence.txt).

Чтение файлов _sequence.txt (FASTQ)

_sequence.txt Файлы - это файлы в формате FASTQ, которые содержат считывание последовательности и их оценки качества после обрезки качества и фильтрации. Вы можете использовать fastqinfo для отображения сводки содержимого _sequence.txt файл и fastqread для чтения содержимого файла. Выходные данные, reads, является массивом ячеек структур, содержащих Header, Sequence и Quality поля.

filename = 'ilmnsolexa_sequence.txt';
info = fastqinfo(filename)
reads = fastqread(filename)

info = 

  struct with fields:

           Filename: 'ilmnsolexa_sequence.txt'
           FilePath: '/tmp/BR2021ad_1584584_202060/publish_examples2/tp64dafef8/ex25447385'
        FileModDate: '06-May-2009 16:02:48'
           FileSize: 30124
    NumberOfEntries: 260


reads = 

  1x260 struct array with fields:

    Header
    Sequence
    Quality

Поскольку на плитку приходится по одному файлу последовательности, нередко в целом имеется коллекция из более чем 1000 файлов. Можно прочесть всю коллекцию файлов, связанных с данным анализом, путем объединения _sequence.txt файлов в один файл. Однако, поскольку эта операция обычно создает большой файл, который требует большого объема памяти для хранения и обработки, рекомендуется считывать содержимое порциями с помощью blockread вариант fastqread функция. Например, можно прочитать первые M последовательностей, или последние M последовательностей, или любые M последовательностей в файле.

M = 150;
N = info.NumberOfEntries;
readsFirst = fastqread(filename, 'blockread', [1 M])
readsLast = fastqread(filename, 'blockread', [N-M+1, N])

readsFirst = 

  1x150 struct array with fields:

    Header
    Sequence
    Quality


readsLast = 

  1x150 struct array with fields:

    Header
    Sequence
    Quality

Съемка распределения длины при считывании последовательности

После загрузки информации о последовательности в рабочую область можно определить количество и длину считываний последовательности и построить график их распределения следующим образом:

seqs = {reads.Sequence};
readsLen = cellfun(@length, seqs);

figure(); hist(readsLen);
xlabel('Number of bases'); ylabel('Number of sequence reads');
title('Length distribution of sequence reads')

Как и ожидалось, в этом примере все операции считывания последовательности имеют длину 36 п.н.

Съемка базового состава считываний последовательности

Вы также можете исследовать нуклеотидный состав, исследуя количество появлений каждого типа основания в каждой считанной последовательности, как показано ниже:

nt = {'A', 'C', 'G', 'T'};
pos = cell(4,N);

for i = 1:4
	pos(i,:) = strfind(seqs, nt{i});
end

count = zeros(4,N);
for i = 1:4
	count(i,:) = cellfun(@length, pos(i,:));
end

%=== plot nucleotide distribution
figure();
subplot(2,2,1); hist(count(1,:)); title('A'); ylabel('Number of sequence reads');
subplot(2,2,2); hist(count(2,:)); title('C');
subplot(2,2,3); hist(count(3,:)); title('G'); xlabel('Occurrences'); ylabel('Number of sequence reads');
subplot(2,2,4); hist(count(4,:)); title('T'); xlabel('Occurrences');

figure(); hist(count');
xlabel('Occurrences');
ylabel('Number of sequence reads');
legend('A', 'C', 'G', 'T');
title('Base distribution by nucleotide type');

Обзор распределения показателей качества

Каждая последовательность считывается в _sequence.txt файл связан с результатом. Балл определяется как Λ = -10 * log10 (p/( 1-p)), где p - ошибка вероятности основания. Можно проверить показатели качества, связанные с базовыми вызовами, преобразовав формат ASCII в числовое представление, а затем построив график их распределения, как показано ниже.

sq = {reads.Quality}; % in ASCII format
SQ = cellfun(@(x) double(x)-64, {reads.Quality}, 'UniformOutput', false); % in integer format

%=== average, median and standard deviation
avgSQ = cellfun(@mean, SQ);
medSQ = cellfun(@median, SQ);
stdSQ = cellfun(@std, SQ);

%=== plot distribution of median and average quality
figure();
subplot(1,2,1); hist(medSQ);
xlabel('Median Score SQ'); ylabel('Number of sequence reads');
subplot(1,2,2); boxplot(avgSQ); ylabel('Average Score SQ');

Преобразование показателей качества между стандартами

Оценки качества, обнаруженные в файлах Solexa/Illumina, асимптотичны, но не идентичны оценкам качества, используемым в стандарте Sanger (Phred-подобные оценки, Q). Q определяется как -10 * log10 (p), где p - вероятность ошибки основания. Например, если показатель качества основы равен Q = 20, то p = 10 ^ (-20/10) = 0,01. Это означает, что на каждые 100 базовых вызовов приходится один неправильный базовый вызов с оценкой 20.

В то время как оценки качества Phred являются положительными целыми числами, оценки качества Solexa/Illumina могут быть отрицательными. Мы можем преобразовать оценки качества Solexa в оценки качества Phred, используя следующий код:

%=== convert from Solexa to Sanger standard
Q = cellfun(@(x) floor(.499 + 10 * log10(1+ 10 .^ (x/10))), SQ, ...
	'UniformOutput', false); % in integer format
q = cellfun(@(x) char(x+33), Q, 'UniformOutput', false); % in ascii format

sanger = q(1:3)'
solexa = sq(1:3)'

sanger =

  3x1 cell array

    {'>>>>>>>>>>>>:>:>>>>>>>>>>>>7&*7.1-%4'}
    {'>>>>>>>>>>>>:>>>>>>>>>:17>5><1;1+&&,'}
    {'>>>>:>>>>>7>5>>>>>5>>>>>7>5.+'69'(-%'}


solexa =

  3x1 cell array

    {']]]]]]]]]]]]Y]Y]]]]]]]]]]]]VCHVMPLAS'}
    {']]]]]]]]]]]]Y]]]]]]]]]YPV]T][PZPICCK'}
    {']]]]Y]]]]]V]T]]]]]T]]]]]V]TMJEUXEFLA'}

Фильтрация и маскирование по показателям качества

Фильтрация чистоты сигнала уже применена к последовательностям в _sequence.txt файлы. Можно выполнить дополнительную фильтрацию, например, рассматривая только те чтения последовательности, базы которых имеют все оценки качества выше определенного порога:

%=== find sequence reads whose bases all have quality above threshold
len = 36;
qt = 10; % minimum quality threshold
a = cellfun(@(x) x > qt, SQ, 'UniformOutput', false);
b = cellfun(@sum, a);
c1 = find(b == len);
n1= numel(c1); % number of sequence reads passing the filter

disp([num2str(n1) ' sequence reads have all bases above threshold ' num2str(qt)]);

30 sequence reads have all bases above threshold 10

Кроме того, можно учитывать только те считывания последовательности, которые имеют меньше заданного числа оснований с показателями качества ниже порогового значения:

%=== find sequence reads having less than M bases with quality below threshold
M = 5; % max number of bases with poor quality
a = cellfun(@(x) x <= qt, SQ, 'UniformOutput', false);
b = cellfun(@sum, a);
c2 = find(b <= M);
n2 = numel(c2); % number of sequence reads passing the filter

disp([num2str(n2) ' sequence reads have less than ' num2str(M) ' bases below threshold ' num2str(qt)]);

235 sequence reads have less than 5 bases below threshold 10

Наконец, маску нижнего регистра можно применить к основам с показателями качества ниже порогового значения:

seq = reads(1).Sequence
mseq = seq;
qt2 = 20;  % quality threshold
mask = SQ{1} < qt;
mseq(mask) = lower(seq(mask))

seq =

    'GGACTTTGTAGGATACCCTCGCTTTCCTTCTCCTGT'


mseq =

    'GGACTTTGTAGGATACCCTCGCTTTCCTtcTCCTgT'

Сводка вхождений для чтения

Для суммирования вхождений чтения можно определить количество уникальных последовательностей чтения и их распределение по набору данных. Можно также идентифицировать те считывания последовательности, которые происходят несколько раз, часто потому, что они соответствуют адаптерам или праймерам, используемым в процессе секвенирования.

%=== determine read frequency
[uReads,~,n] = unique({reads.Sequence});
numUnique = numel(uReads)
readFreq = accumarray(n(:),1);
figure(); hist(readFreq, unique(readFreq));
xlabel('Occurrences'); ylabel('Number of sequence reads');
title('Read occurrences');

%=== identify multiply-occurring sequence reads
d = readFreq > 1;
dupReads = uReads(d)'
dupFreq  = readFreq(d)'

numUnique =

   250


dupReads =

  9x1 cell array

    {'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'}
    {'GATTTTATTGGTATCAGGGTTAATCGTGCCAAGAAA'}
    {'GCATGGGTGATGCTGGTATTAAATCTGCCATTCAAG'}
    {'GGGATGAACATAATAAGCAATGACGGCAGCAATAAA'}
    {'GGGGGAGCACATTGTAGCATTGTGCCAATTCATCCA'}
    {'GGTTATTAAAGAGATTATTTGTCTCCAGCCACTTAA'}
    {'GTTCTCACTTCTGTTACTCCAGCTTCTTCGGCACCT'}
    {'GTTGCTGCCATCTCAAAAACATTTGGACTGCTCCGC'}
    {'GTTGGTTTCTATGTGGCTAAATACGTTAACAAAAAG'}


dupFreq =

     2     2     2     2     2     2     3     2     2

Идентификация артефактов гомополимеров

Секвенирование Illumina/Solexa может создавать ложные поли на краях плитки. Чтобы идентифицировать эти артефакты, нужно идентифицировать гомополимеры, то есть считывает последовательность, составленную только из одного типа нуклеотидов. В рассматриваемом наборе данных имеется два гомополимера, оба из которых являются полиА.

%=== find homopolymers
pc = (count ./ len) * 100;
[homopolType,homopolIndex] = find(pc == 100);

homopolIndex
homopol = {reads(homopolIndex).Sequence}'

homopolIndex =

   251
   257


homopol =

  2x1 cell array

    {'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'}
    {'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'}

Аналогично можно идентифицировать считывания последовательностей, которые близки к гомополимерам, то есть считывания последовательностей, которые состоят почти исключительно из одного нуклеотидного типа.

%=== find near-homopolymers
[nearhomopolType, nearhomopolIndex] = find(pc < 100 & pc > 85); % more than 85% same base
nearhomopolIndex
nearHomopol = {reads(nearhomopolIndex).Sequence}'

nearhomopolIndex =

     4
   243


nearHomopol =

  2x1 cell array

    {'AAAAACATAAAAAAAAAAATAAAAAAACAAAAAAAA'}
    {'AAAAAAATAAAAAAAAAAATAAAAAAAAATTAAAAA'}

Запись данных в формат FASTQ

После обработки и анализа данных может оказаться удобным сохранить подмножество последовательностей в отдельном файле FASTQ для дальнейшего рассмотрения. Для этого можно использовать fastqwrite функция.