Этот пример показывает, как считать и выполнить основные операции с данными, произведенными Illumina/Solexa Genome Analyzer®.
Во время анализа, запущенного с программным обеспечением Genome Analyzer Pipeline, производятся несколько промежуточных файлов. В этом примере вы изучите, как считать и управлять информацией, содержавшейся в файлах последовательности (_sequence.txt).
Файлы _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/BR2019ad_1035872_198992/publish_examples4/tpe1748eff/ex25447385'
FileModDate: '06-May-2009 16:02:48'
FileSize: 30124
NumberOfEntries: 260
reads =
1x260 struct array with fields:
Header
Sequence
Quality
Поскольку существует один файл последовательности на мозаику, весьма распространено иметь набор более чем 1 000 файлов всего. Можно считать целый набор файлов, сопоставленных с данным аналитическим шансом путем конкатенации файлов _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 BP долго.
Можно также исследовать состав нуклеотида путем рассмотрения количества случаев каждого базового типа в каждом чтении последовательности, как показано ниже:
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 сопоставлено со счетом. Счет задан как SQ =-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), =.01. Это означает, что существует один неправильный базовый вызов каждые 100 основных вызовов со счетом of20.
В то время как качественные очки 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 для будущего фактора. С этой целью можно использовать функцию fastqwrite.