В этом примере показано, как считать и выполнить основные операции с данными, произведенными 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/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
Поскольку существует один файл последовательности на мозаику, весьма распространено иметь набор более чем 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
функция.