На этой странице выложен отчет по практикуму 15.
Задание 1. Подготовка чтений программой trimmomatic.
На странице проекта http://www.ebi.ac.uk/ena/data/view/SRR4240381 скачан файл fastq в формате .gz.
Скачанный файл перенесен в рабочую директорию /nfs/srv/databases/ngs/tsvetcovroman. Там он был распакован программой gunzip в командной строке с помощью PuTTy.
При этом использована команда:
gunzip SRR4240380.fastq.gz
В результате получен файл SRR4240380.fastq.
Далее из файла были удалены остатки адаптеров. В качестве адаптеров использованы адаптеры из файлов в директории /P/y15/term3/block4/adapters.
Эти файлы были объединены в файл ad.fasta в директории /nfs/srv/databases/ngs/tsvetcovroman.
Для удаления остатков адаптеров использована команда:
java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240380.fastq SRR4240380_out.fastq ILLUMINACLIP:ad.fasta:2:7:7
Вывод программы:
TrimmomaticSE: Started with arguments: -phred33 SRR4240380.fastq SRR4240380_out.fastq ILLUMINACLIP:ad.fasta:2:7:7
Automatically using 8 threads
Using PrefixPair: 'AGATGTGTATAAGAGACAG' and 'AGATGTGTATAAGAGACAG'
Using PrefixPair: 'TACACTCTTTCCCTACACGACGCTCTTCCGATCT' and 'GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT'
Using Long Clipping Sequence: 'GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'TTTTTTTTTTAATGATACGGCGACCACCGAGATCTACAC'
Using Long Clipping Sequence: 'TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'TTTTTTTTTTCAAGCAGAAGACGGCATACGA'
Using Long Clipping Sequence: 'CTGTCTCTTATACACATCTGACGCTGCCGACGA'
Using Long Clipping Sequence: 'AGATCGGAAGAGCTCGTATGCCGTCTTCTGCTTG'
Using Long Clipping Sequence: 'AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC'
Using Long Clipping Sequence: 'TACACTCTTTCCCTACACGACGCTCTTCCGATCT'
Using Long Clipping Sequence: 'GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT'
Using Long Clipping Sequence: 'AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA'
Using Long Clipping Sequence: 'AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT'
Using Long Clipping Sequence: 'AGATCGGAAGAGCGGTTCAGCAGGAATGCCGAG'
Skipping duplicate Clipping Sequence: 'AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC'
Using Long Clipping Sequence: 'AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT'
Using Long Clipping Sequence: 'AGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGTCTTCTGCTTG'
Using Long Clipping Sequence: 'CAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATCT'
Skipping duplicate Clipping Sequence: 'AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA'
Using Long Clipping Sequence: 'CTGTCTCTTATACACATCTCCGAGCCCACGAGAC'
Using Long Clipping Sequence: 'AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT'
ILLUMINACLIP: Using 2 prefix pairs, 17 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Input Reads: 5217318 Surviving: 5119144 (98,12%) Dropped: 98174 (1,88%)
TrimmomaticSE: Completed successfully
В результате был получе файл SRR4240380_out.fastq.
До очиски чтений и удаления плохих букв с концов чтений , размер файла
SRR4240380.fastq составил 538429 (526M ). Для получения этих результатов использованы команды ls -s(выдает результат 538429) и du -h(выдает результат 526M).
Далее были удалены плохие буквы с концов чтений (как в практикуме 13), и оставлены только чтения длиной не менее 30.
Для этого использовалась команда:
java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240380_out.fastq SRR.fastq TRAILING:20 MINLEN:30
Выдача команды:
TrimmomaticSE: Started with arguments: -phred33 SRR4240380_out.fastq SRR.fastq TRAILING:20 MINLEN:30
Automatically using 8 threads
Input Reads: 5119144 Surviving: 4879709 (95,32%) Dropped: 239435 (4,68%)
TrimmomaticSE: Completed successfully
В результате получен файл SRR.fastq. Размер этого файл меньше размера исходного файла и составляет 513137(данные получены с помощью команды
ls -s). Размер файла сократился на 4,7%.
Затем с помощью команды fastqc SRR4240380_out.fastq
был проанализирован файл SRR4240380_out.fastq.
Результаты:SRR4240380_out_fastqc.html. В нем 5119144 чтений.
Затем с помощью команды fastqc SRR.fastq
был проанализирован файл SRR.fastq.
Результаты:SRR_fastqc.html. В нем 4981466 чтений.
То есть в результате работы программы было удалено 235852 ридов, что составляет 4,52% от общего количества.
Далее была запущена программа velveth, чтобы подготовить k-меры длины k=29 (максимально возможной при нашей длине чтений). Файлы записывались в директорию dir.
k-меры длины k=29 были подготовлены с помощью команды:
velveth dir 29 -short -fastq SRR.fastq
Выдача программы:
[0.000000] Reading FastQ file SRR.fastq
[18.396233] 4981466 reads found.
[18.396282] Done
[18.491252] Reading read set file dir/Sequences;
[21.265722] 4981466 sequences found
[33.830727] Done
[33.830789] 4981466 sequences in total.
[33.830962] Writing into roadmap file dir/Roadmaps...
[38.900422] Inputting sequences...
[38.944590] Inputting sequence 0 / 4981466
[48.647606] Inputting sequence 2000000 / 4981466
[56.191919] Inputting sequence 4000000 / 4981466
[61.957199] Inputting sequence 1000000 / 4981466
[79.406068] Inputting sequence 3000000 / 4981466
[94.728699] === Sequences loaded in 56.664041 s
[94.728932] Done inputting sequences
[94.728938] Destroying splay table
[94.963624] Splay table destroyed
Сборка на основе k-меров производилась при помощи программы velvetg, которая строит и анализирует граф де-Брёйна.
В нём k-меры являются вершинами, а ориентированные рёбра соединяют так называемые префикс- и суффикс-последовательности.
При первом запуске авторы рекомендуют указывать параметр порогового покрытия так: -cov_cutoff auto.
Значение auto рассчитывается как половина от длины взвешенной медианы по глубине покрытия контига. Итоговая команда:
velvetg dir -cov_cutoff auto
Выдача программы:
[0.000001] Reading roadmap file dir/Roadmaps
[19.017692] 4981466 roadmaps read
[19.041832] Creating insertion markers
[26.771227] Ordering insertion markers
[30.191986] Counting preNodes
[32.839218] 1023268 preNodes counted, creating them now
[38.237750] Sequence 1000000 / 4981466
[43.225072] Sequence 2000000 / 4981466
[47.382583] Sequence 3000000 / 4981466
[52.255024] Sequence 4000000 / 4981466
[56.573892] Adjusting marker info...
[59.455779] Connecting preNodes
[60.572649] Connecting 2000000 / 4981466
[61.770981] Connecting 4000000 / 4981466
[62.785777] Connecting 1000000 / 4981466
[63.800992] Connecting 3000000 / 4981466
[65.136966] Cleaning up memory
[65.182918] Done creating preGraph
[65.182963] Concatenation...
[66.685091] Renumbering preNodes
[66.685136] Initial preNode count 1023268
[66.839233] Destroyed 746274 preNodes
[66.839276] Concatenation over!
[66.839288] Clipping short tips off preGraph
[67.341662] Concatenation...
[67.898046] Renumbering preNodes
[67.898359] Initial preNode count 276994
[67.904486] Destroyed 268912 preNodes
[67.916921] Concatenation over!
[67.916975] 193620 tips cut off
[67.917024] 8082 nodes left
[67.917254] Writing into pregraph file dir/PreGraph...
[68.160901] Reading read set file dir/Sequences;
[70.498173] 4981466 sequences found
[82.597148] Done
[86.509429] Reading pre-graph file dir/PreGraph
[87.318856] Graph has 8082 nodes and 4981466 sequences
[87.433243] Scanning pre-graph file dir/PreGraph for k-mers
[87.622596] 697680 kmers found
[87.809551] Sorting kmer occurence table ...
[88.422370] Sorting done.
[88.422591] Computing acceleration table...
[88.543097] Computing offsets...
[88.571004] Ghost Threading through reads 2000000 / 4981466
[88.578979] Ghost Threading through reads 3000000 / 4981466
[88.584186] Ghost Threading through reads 1000000 / 4981466
[88.585342] Ghost Threading through reads 4000000 / 4981466
[88.598005] Ghost Threading through reads 0 / 4981466
[88.625380] === Ghost-Threaded in 0.073996 s
[88.644901] Threading through reads 0 / 4981466
[90.408260] Threading through reads 2000000 / 4981466
[92.556389] Threading through reads 4000000 / 4981466
[94.557830] Threading through reads 1000000 / 4981466
[96.308401] Threading through reads 3000000 / 4981466
[98.403365] === Threaded in 9.777908 s
[98.498046] Correcting graph with cutoff 0.200000
[98.499313] Determining eligible starting points
[98.507633] Done listing starting nodes
[98.516776] Initializing todo lists
[98.517859] Done with initilization
[98.517880] Activating arc lookup table
[98.518930] Done activating arc lookup table
[98.571505] Concatenation...
[98.571906] Renumbering nodes
[98.571920] Initial node count 8082
[98.572162] Removed 386 null nodes
[98.572176] Concatenation over!
[98.572184] Clipping short tips off graph, drastic
[98.574405] Concatenation...
[98.598886] Renumbering nodes
[98.598911] Initial node count 7696
[98.599023] Removed 6809 null nodes
[98.599041] Concatenation over!
[98.599055] 887 nodes left
[98.599219] Writing into graph file dir/Graph...
[98.865882] Measuring median coverage depth...
[98.866235] Median coverage depth = 34.814990
[98.866952] Removing contigs with coverage < 17.407495...
[98.870778] Concatenation...
[98.874402] Renumbering nodes
[98.888899] Initial node count 887
[98.888945] Removed 745 null nodes
[98.888966] Concatenation over!
[98.888992] Concatenation...
[98.889017] Renumbering nodes
[98.889031] Initial node count 142
[98.889047] Removed 0 null nodes
[98.889060] Concatenation over!
[98.889074] Clipping short tips off graph, drastic
[98.889108] Concatenation...
[98.889177] Renumbering nodes
[98.889196] Initial node count 142
[98.889224] Removed 51 null nodes
[98.889238] Concatenation over!
[98.889256] 91 nodes left
[98.889267] WARNING: NO EXPECTED COVERAGE PROVIDED
[98.889286] Velvet will be unable to resolve any repeats
[98.889297] See manual for instructions on how to set the expected coverage parameter
[98.889313] Concatenation...
[98.889332] Renumbering nodes
[98.889351] Initial node count 91
[98.889363] Removed 0 null nodes
[98.889380] Concatenation over!
[98.889396] Removing reference contigs with coverage < 17.407495...
[98.889426] Concatenation...
[98.889440] Renumbering nodes
[98.889459] Initial node count 91
[98.889478] Removed 0 null nodes
[98.889494] Concatenation over!
[98.965326] Writing contigs into dir/contigs.fa...
[99.147753] Writing into stats file dir/stats.txt...
[99.148301] Writing into graph file dir/LastGraph...
[99.411193] Estimated Coverage cutoff = 17.407495
Final graph has 91 nodes and n50 of 26659, max 64995, total 655256, using 0/4981466 reads
Результаты работы velvetg:
изначально 4981466 вершин, в финальной версии графа — 91;
•N50 = 26659;
•файлы с графами: PreGraph, LastGraph, Graph;
•stats.txt — информация о вершинах графа: 1.длина вершины в k-мерах (для расчёта длины в нуклеотидах нужно просто добавить k–1);
2.in — число входящих рёбер с 5'-конца;
3.out — число выходящих рёбер с 3'-конца;
4.колонки *_cov и *_Ocov несут информацию о покрытии, различаются между собой методом вычислений (*_Ocov — более строгие расчёты).
•contigs.fa — [fasta]-файл с контигами — последовательностями длины не менее 2k (мы не устанавливали специальное ограничение
на минимально возможную длину контига, но такая возможность, вообще говоря, существует).
•анализ контигов удобнее проводить по stats.txt: конкретные нуклеотидные последовательности нас пока не особо интересуют.
При этом надо учитывать, что не все вершины стали контигами (всего их 48 штук, а вершин, как было сказано ранее, 91).
Нужно сперва отсечь все те, длина которых меньше 29, ибо они не прописываются в contigs.fa и не являются "полноценными контигами".
•данные по контигам: среднее значение длины контигов (в k-мерах): 7200,615;
медианное значение длины: 224;
файл с расчетами
самые длинные контиги: