Введение

Во время анализа, запущенного с программным обеспечением 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/BR2021bd_1724986_5107/publish_examples2/tp6f4f0efe/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

Если вы обработали и анализировали свои данные, может быть удобно сохранить подмножество последовательностей в отдельном файле FASTQ для будущего фактора. С этой целью можно использовать fastqwrite функция.