Сначала чтения были скачаны с помощью следующей команды:
$ wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR424/006/SRR4240356/SRR4240356.fastq.gz
Чтобы было удобнее удалять адаптеры, их можно собрать в один файл:
$ cat /mnt/scratch/NGS/adapters/* >> adapters.fasta
Чтобы удалить адаптеры, была применим следующую команду:
$ java -jar /usr/share/java/trimmomatic.jar SE -threads 10 -phred33 -trimlog trimlog.txt SRR4240356.fastq.gz trim.fastq.gz ILLUMINACLIP:adapters.fasta:2:7:7
Таким образом было удалено 2.04% чтений.
Далее нужно удалить с правых концов чтений такие нуклеотиды, которые имеют качество меньше 20, а также удалить чтения длиной меньше 32. Снова применим trimmomatic:
$ java -jar /usr/share/java/trimmomatic.jar SE -threads 10 -phred33 -trimlog trimlog2.txt trim.fastq.gz trim2.fastq.gz TRAILING:20 MINLEN:32
После этого было удалено 4.15% чтений. Нетрудно посчитать, что после обеих операций было оставлено 93.89% чтений. Исходный файл весил 174.26 мб, а самый последний получившийся весит 162.02 мб.
Поулчение k-меров и сборка генома
Сначала нужно получить k-меры. Для этого применим следующую команду:
$ velveth assemble 31 -fastq.gz -short trim2.fastq.gz &> velveth.log
В директории assemble было создано три файла: Log, Roadmaps, Sequences. Теперь можно приступить к сборке. Для этого применим следующую команду:
$ velvetg assemble &> velvetg.log
Полученная сборка имеет N50 равное 65554 (т.е. 65584 нуклеотидов). Чтобы узнать, какие из контигов самые длинные, можно применить такую команду:
Сильнее всех выделяется контиг с ID 64, который имеет покрытие 266951.0, что во много раз больше покрытия любого другого контига, но он имеет длину 1 (т.е. 31 нуклеотид). Из неединичных контигов лучше всего покрыт контиг с ID 27, и его покрытие составляет 458.43. Что касается слабопокрытых участков, то есть несколько контигов с покрытием ровно 1.0, все они имеют длину примерно 1 и перечислены в выдаче выше.
Чтобы достать последовательности наибольших контигов, применим seqret:
$ seqret assemble/contigs.fa:'*length_111962*' contig8.fasta
Read and write (return) sequences
$ seqret assemble/contigs.fa:'*length_107488*' contig6.fasta
Read and write (return) sequences
$ seqret assemble/contigs.fa:'*length_80939*' contig10.fasta
Read and write (return) sequences
Анализ полученных контигов
Чтобы понять, как полученные контиги лежат на хромосоме Buchnera aphidicola (AC: CP009253) воспользуемся megablast для двух последовательностей. Выдача BLAST для каждого из контигов состоит из большого количества выравниваний (от 11 до 18), и ложатся все контиги с разрывами. Тем не менее, выдачу можно отсортировать по координатам в последовательности хромосомы, поэтому левой координатой контига на хромосоме назовём левую координату самого первого выравнивания, а правой — правую координату самого последнего. К сожалению, при скачивании выдачи сортировка сбивается.
Все контиги ложатся на хромосому без перестановок, т.е. чем правее нуклеотид, тем правее он ляжет на хромосоме, но, опять же, имеются разрывы. Причём все эти разрывы — это не "делеции", а именно участки, не похожие друг на друга в контиге и хромосоме, примерно равной длины. Контиг 8 лежит на хромосоме по координатам 451729 — 555905. Суммарный вес выравниваний составляет 50932. Суммарное Identity выравнивания контига с хромосомой составляет 81.46%. Покрытие равно 75%. Имеется 5 крупных разрывов разной длины, но самый крупный имеет длину примерно 10 000 нуклеотидов. Карту локального сходства можно видеть на Рисунке 1. С выдачей BLAST можно ознакомиться по ссылке. Рисунок 1. Карта локального сходства контига 8 и хромосомы Buchnera aphidicola.Контиг 6 лежит на хромосоме по координатам 220869 — 223720. Суммарный вес выравниваний составляет 43281. Суммарное Identity выравнивания контига с хромосомой составляет 78.76%. Покрытие равно 74%. Имеется 9 крупных разрывов, самый длинный из которых составляет примерно 7000 нуклеотидов. Карту локального сходства можно видеть ниже на Рисунке 2. С выдачей BLAST можно ознакомиться по ссылке. Рисунок 2. Карта локального сходства контига 6 и хромосомы Buchnera aphidicola.Контиг 10 лежит на хромосоме по координатам 126623 — 195400. Суммарный вес выравниваний равен 30801. Суммарное Identity выравнивания контига с хромосомой составляет 81.46%. Покрытие равно 74.88%. Имеется 5 разрывов длиной не больше 7000 нуклеотидов. Карту локального сходства можно видеть на Рисунке 3. Интересно, что полученный контиг обратно комплементарен последовательности хромосомы. С выдачей BLAST можно ознакомиться по ссылке. Рисунок 3. Карта локального сходства контига 10 и хромосомы Buchnera aphidicola.
