Введение

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

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

The _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 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 может давать ложный polyA на ребрах плитки. Чтобы идентифицировать эти программные продукты, необходимо идентифицировать гомополимеры, то есть последовательности, состоящие только из одного типа нуклеотидов. В рассматриваемом наборе данных существует два гомополимера, оба из которых являются полиА.

%=== 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 функция.