Class04_NGS
Alfredo Cárdenas-Rivera
9/10/2020
SECUENCIACIÓN DE NUEVA GENERACIÓN
Importando secuencias con al función readFast() del paquete “ShortRead”. Este paquete permite la importancpón de secuencias en diferentes formatos. Primero llamamos a la librería de “shortRead”, que demora algunos minutos en instalarse
library(ShortRead)## Loading required package: BiocGenerics## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which.max, which.min## Loading required package: BiocParallel## Loading required package: Biostrings## Loading required package: S4Vectors## Loading required package: stats4##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:base':
##
## expand.grid## Loading required package: IRanges##
## Attaching package: 'IRanges'## The following object is masked from 'package:grDevices':
##
## windows## Loading required package: XVector##
## Attaching package: 'Biostrings'## The following object is masked from 'package:base':
##
## strsplit## Loading required package: Rsamtools## Loading required package: GenomeInfoDb## Loading required package: GenomicRanges## Loading required package: GenomicAlignments## Loading required package: SummarizedExperiment## Loading required package: MatrixGenerics## Loading required package: matrixStats##
## Attaching package: 'MatrixGenerics'## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attaching package: 'Biobase'## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians1. IMPORTACIÓN DE BASES DE DATOS (FASTQ)
Una vez instalado el paquete, enlistamos todos los archivos a analizar
pwd <- paste0(dirname(getwd()), "/data/clase_4/brca/")
print(pwd)## [1] "C:/Bioinformatica_2020_2_BI0035/data/clase_4/brca/"fastaqFiles <- list.files(pwd)
print(fastaqFiles)## [1] "brca.example.illumina.0.1.fastq"
## [2] "brca.example.illumina.0.1.fastq_filter.fastq"
## [3] "brca.example.illumina.0.2.fastq"
## [4] "brca.example.illumina.bowtie2.bam"
## [5] "brca.example.illumina.hg18.vcf"Como vemos, la carpeta brca está constituida por 4 archivos:
- 2 archivos .fastq (brca.example.illumina.0.1.fastq y brca.example.illumina.0.2.fastq)
- 1 archivo .bam (brca.example.illumina.bowtie2.bam)
- 1 archivo .vcf (brca.example.illumina.hg18.vcf)
Para importar todas los archivos a un solo objeto lista en R, podemos usar la función readFastq() del paquete “ShortRead”. Esta función tiene 3 parámetros:
- dirPath = donde se define la carpeta donde se encuentran los archivos
- pattern = donde se define el patron de los archivos a ser descargados
- full.names = donde se define si se desea conservar el nombre completo original de los archivos importados.
fq <- ShortRead::readFastq(dirPath = pwd, pattern = "*.fastq",full.names = T)
fq## class: ShortReadQ
## length: 93172 reads; width: 81..101 cyclesVemos que el objeto creado “fq” es de la clase “ShortReadQ” y que tiene una longitud de 63426 lecturas en 101 ciclos.
2. EXPLORACIÓN DE LOS ARCHIVOS FASTQ IMPORTADOS
Para explorar los Ids de los archivos .fastq importados podemos usar la función ShortRead::id():
ShortRead::id(fq)## BStringSet object of length 93172:
## width seq
## [1] 55 Frag_1 chr17 (Strand + Offset 106709--107175) 467M 101M
## [2] 51 Frag_2 chr17 (Strand + Offset 5449--5912) 464M 101M
## [3] 53 Frag_3 chr17 (Strand + Offset 82312--82850) 539M 101M
## [4] 51 Frag_4 chr17 (Strand - Offset 1052--1590) 539M 101M
## [5] 51 Frag_5 chr17 (Strand + Offset 3003--3748) 746M 101M
## ... ... ...
## [93168] 57 Frag_15854 chr17 (Strand - Offset 45380--46072) 693M 101M
## [93169] 57 Frag_15855 chr17 (Strand - Offset 17933--18549) 617M 101M
## [93170] 57 Frag_15856 chr17 (Strand + Offset 65810--66418) 609M 101M
## [93171] 59 Frag_15857 chr17 (Strand - Offset 116176--116906) 731M 101M
## [93172] 59 Frag_15858 chr17 (Strand - Offset 116176--116906) 731M 101MPara explorar las secuencias podemos usar la función ShortRead::sread():
ShortRead::sread(fq)## DNAStringSet object of length 93172:
## width seq
## [1] 101 AGGCAGGTCTCAAACTCCTGACCTCAGGTGAT...GCATGAGCCACCATGTCCGGCAAGTTTCTTT
## [2] 101 ACTCCTGACAAGTGATCCACCTGCCTCGGCCT...CAGCCTCCAGCCCATCATTTCTTGATGATTT
## [3] 101 TTTGCCATAAAGTGCCTGCCCTCTAGCCTCTA...GTACGAAGGTCAGAATCGCTACCTATTGTCC
## [4] 101 TTGCAGTGAGCCAAGATCATACCACGGCACTC...AAAAAAAAAAGGAAAATGAAACTAGAAGAGA
## [5] 101 CTCAGCTCACTGCAACCTCCACCTCCCGGGTT...ATTAAAGGTGCCTGCCACCATGCCAGGCTAA
## ... ... ...
## [93168] 101 TACTAGCAGCTAAAAGGGTATTTAAGATTTAG...TCTGAACAGCTTCATTTTAAGAATGTGCCTA
## [93169] 101 TTCAAGGACAATTATTTCTTGTTATTTAGGCG...TTAAACGATAGTACTCACCAACCCTTGCCAG
## [93170] 101 TGATGCCTTCTTGATCCATTATACTAGAGTAT...ATTCTGAGTTTTAAAGGTACTTATGTACAGT
## [93171] 101 ACGCTCCTCACTTCCCAGACGGGGTGGCTGCC...CCTGGTAGAGGGTCTCGTCACTTCTCAGACG
## [93172] 101 ACGCTCCTCACNNCCCAGACGGGGTGNNNGCC...CCTGGTAGAGGGTCTCGTNNNNNNTCAGACGY para ver la calidad de las lecturas podemos usar la función Biostrings::quality():
Biostrings::quality(fq)## class: FastqQuality
## quality:
## BStringSet object of length 93172:
## width seq
## [1] 101 CGED?FG?GGFFGFGFBFFGEGGEGFEGFGCG...GGG<GEGGF@GGFGGEEG4G5DFFF@F8GGD
## [2] 101 DCGCGGGGFGGFEGGGC=FFFGEGGGGDFG=F...FGGBEGFFAGEBFGGFGGGGCDGFD;FGG>G
## [3] 101 GGGG7F@DGGGFDGGGGGGAG7FGGFGGGFGF...GFCGGFEGGGFG>GGGFD)DFAFGDGCGFG5
## [4] 101 GGEGG?GGGFGEGGG@*F@FFFEGGF>G?GGF...FGEFGCGGFGGGFCCFGGFFF6?E=;6B@FC
## [5] 101 GGFGFEFCFEFG<GGG>82F46GGGEEGC8GG...G0FC(F%FGGC3EFFDAD'9=;16BBGE@<3
## ... ... ...
## [93168] 101 A=GGG@GBGGG5%DF?GDEA?DE?GGEEEG>F...FDG>5,C,?2F9>BDB98>?F>*-=E:,310
## [93169] 101 ;A<7C2>>8A?6A,@8(C8*;>=3D;9,4,);...7:5$;;.0(2>=&-<3*&7+6(1(&1-&1.2
## [93170] 101 FGGGGGFGGFGGGGGGGGGGEGGFGBGGBDGG...F8GGDEEFG*@GG%FGFEDFGDG@CFB51=F
## [93171] 101 GFFG5(DBGGF4EC2FFGFGEGAFFD<GEGDA...<B#BBE8E:BED4==A&E<?;ED>=:1;:25
## [93172] 101 GFFG5(DBGGF4EC2FFGFGEGAFFD<GEGDA...<B#BBE8E:BED4==A&E<+-&++-((&&-+Tambien podemos obtener una visión rápida de las secuencias con la función Biostrings::detail()
Biostrings::detail(fq)## class: ShortReadQ
##
## sread:
## DNAStringSet object of length 93172:
## width seq
## [1] 101 AGGCAGGTCTCAAACTCCTGACCTCAGGTGAT...GCATGAGCCACCATGTCCGGCAAGTTTCTTT
## [2] 101 ACTCCTGACAAGTGATCCACCTGCCTCGGCCT...CAGCCTCCAGCCCATCATTTCTTGATGATTT
## [3] 101 TTTGCCATAAAGTGCCTGCCCTCTAGCCTCTA...GTACGAAGGTCAGAATCGCTACCTATTGTCC
## [4] 101 TTGCAGTGAGCCAAGATCATACCACGGCACTC...AAAAAAAAAAGGAAAATGAAACTAGAAGAGA
## [5] 101 CTCAGCTCACTGCAACCTCCACCTCCCGGGTT...ATTAAAGGTGCCTGCCACCATGCCAGGCTAA
## ... ... ...
## [93168] 101 TACTAGCAGCTAAAAGGGTATTTAAGATTTAG...TCTGAACAGCTTCATTTTAAGAATGTGCCTA
## [93169] 101 TTCAAGGACAATTATTTCTTGTTATTTAGGCG...TTAAACGATAGTACTCACCAACCCTTGCCAG
## [93170] 101 TGATGCCTTCTTGATCCATTATACTAGAGTAT...ATTCTGAGTTTTAAAGGTACTTATGTACAGT
## [93171] 101 ACGCTCCTCACTTCCCAGACGGGGTGGCTGCC...CCTGGTAGAGGGTCTCGTCACTTCTCAGACG
## [93172] 101 ACGCTCCTCACNNCCCAGACGGGGTGNNNGCC...CCTGGTAGAGGGTCTCGTNNNNNNTCAGACG
##
## id:
## BStringSet object of length 93172:
## width seq
## [1] 55 Frag_1 chr17 (Strand + Offset 106709--107175) 467M 101M
## [2] 51 Frag_2 chr17 (Strand + Offset 5449--5912) 464M 101M
## [3] 53 Frag_3 chr17 (Strand + Offset 82312--82850) 539M 101M
## [4] 51 Frag_4 chr17 (Strand - Offset 1052--1590) 539M 101M
## [5] 51 Frag_5 chr17 (Strand + Offset 3003--3748) 746M 101M
## ... ... ...
## [93168] 57 Frag_15854 chr17 (Strand - Offset 45380--46072) 693M 101M
## [93169] 57 Frag_15855 chr17 (Strand - Offset 17933--18549) 617M 101M
## [93170] 57 Frag_15856 chr17 (Strand + Offset 65810--66418) 609M 101M
## [93171] 59 Frag_15857 chr17 (Strand - Offset 116176--116906) 731M 101M
## [93172] 59 Frag_15858 chr17 (Strand - Offset 116176--116906) 731M 101M
## class: FastqQuality
## quality:
## BStringSet object of length 93172:
## width seq
## [1] 101 CGED?FG?GGFFGFGFBFFGEGGEGFEGFGCG...GGG<GEGGF@GGFGGEEG4G5DFFF@F8GGD
## [2] 101 DCGCGGGGFGGFEGGGC=FFFGEGGGGDFG=F...FGGBEGFFAGEBFGGFGGGGCDGFD;FGG>G
## [3] 101 GGGG7F@DGGGFDGGGGGGAG7FGGFGGGFGF...GFCGGFEGGGFG>GGGFD)DFAFGDGCGFG5
## [4] 101 GGEGG?GGGFGEGGG@*F@FFFEGGF>G?GGF...FGEFGCGGFGGGFCCFGGFFF6?E=;6B@FC
## [5] 101 GGFGFEFCFEFG<GGG>82F46GGGEEGC8GG...G0FC(F%FGGC3EFFDAD'9=;16BBGE@<3
## ... ... ...
## [93168] 101 A=GGG@GBGGG5%DF?GDEA?DE?GGEEEG>F...FDG>5,C,?2F9>BDB98>?F>*-=E:,310
## [93169] 101 ;A<7C2>>8A?6A,@8(C8*;>=3D;9,4,);...7:5$;;.0(2>=&-<3*&7+6(1(&1-&1.2
## [93170] 101 FGGGGGFGGFGGGGGGGGGGEGGFGBGGBDGG...F8GGDEEFG*@GG%FGFEDFGDG@CFB51=F
## [93171] 101 GFFG5(DBGGF4EC2FFGFGEGAFFD<GEGDA...<B#BBE8E:BED4==A&E<?;ED>=:1;:25
## [93172] 101 GFFG5(DBGGF4EC2FFGFGEGAFFD<GEGDA...<B#BBE8E:BED4==A&E<+-&++-((&&-+3. EXPLORACIÓN POR SECUENCIAS INDIVIDUALES
Si desamos ver las sesecuencias de forma individual pordmemo usar la función ShortRead::sread(fq)[i]
ShortRead::sread(fq)[256]## DNAStringSet object of length 1:
## width seq
## [1] 101 CCTCTTGAGATGGGTAGTTTCTATTCTGAAGACT...AATGATGGGCATTTAGAAGGGGATGACCTAGAAPodemos extraer la información de la secuencia convirtiendo el resultado de la función a un caracter o a un vector
as.character(ShortRead::sread(fq)[256])## [1] "CCTCTTGAGATGGGTAGTTTCTATTCTGAAGACTCCCAGAGCAACTGTGCATGTACCACCTATCATCTAATGATGGGCATTTAGAAGGGGATGACCTAGAA"as.vector(ShortRead::sread(fq)[256])## [1] "CCTCTTGAGATGGGTAGTTTCTATTCTGAAGACTCCCAGAGCAACTGTGCATGTACCACCTATCATCTAATGATGGGCATTTAGAAGGGGATGACCTAGAA"Podemos serparar cada nucleótido de la secuencia mediante la función strsplit(), definiendo en el parámetro de separación como ""
strsplit(x = as.character(ShortRead::sread(fq)[256]),"")## [[1]]
## [1] "C" "C" "T" "C" "T" "T" "G" "A" "G" "A" "T" "G" "G" "G" "T" "A" "G" "T"
## [19] "T" "T" "C" "T" "A" "T" "T" "C" "T" "G" "A" "A" "G" "A" "C" "T" "C" "C"
## [37] "C" "A" "G" "A" "G" "C" "A" "A" "C" "T" "G" "T" "G" "C" "A" "T" "G" "T"
## [55] "A" "C" "C" "A" "C" "C" "T" "A" "T" "C" "A" "T" "C" "T" "A" "A" "T" "G"
## [73] "A" "T" "G" "G" "G" "C" "A" "T" "T" "T" "A" "G" "A" "A" "G" "G" "G" "G"
## [91] "A" "T" "G" "A" "C" "C" "T" "A" "G" "A" "A"Sobre este resultado podemos hacer algunos cálculos básicos como una tabla de contingencia o un barpolot
x <- strsplit(x = as.character(ShortRead::sread(fq)[256]),"")
table(x)## x
## A C G T
## 28 21 24 28barplot(table(x))
De la misma forma podemos obetener información de la calidad de la secuenciación usando la función Biostring::quality()
x <- strsplit(x = as.character(Biostrings::quality(fq)[[256]]),"")
x## [[1]]
## [1] "F" "F" "G" "F" "G" "G" "F" "G" "G" "G" "2" "G" "G" "E" "E" "@" "@" "G"
## [19] "F" "G" "G" "G" "G" "(" "G" "D" "4" "E" "D" "?" "G" "D" "G" "G" "G" "G"
## [37] "G" "G" "G" "G" "G" "G" "G" "G" "G" "@" "F" "F" "G" "G" "F" "F" "G" "G"
## [55] "G" "F" "F" "G" "C" "G" "G" "=" "G" "F" "9" "F" "F" "G" "E" "*" "G" "G"
## [73] "E" ":" "G" "F" "G" "F" "B" "G" "G" "G" "@" "E" "G" "6" "G" "G" "F" "F"
## [91] "," "F" "F" ";" "E" "C" "A" "E" "F" ">" "@"table(x)## x
## ( * , : ; ? @ = > 2 4 6 9 A B C D E F G
## 1 1 1 1 1 1 5 1 1 1 1 1 1 1 1 2 3 8 21 48====================================================================================================================================
4. SIMULACIÓN DE DATOS PARA ESTUDIO
Muchas veces la cantidad de información de lecturas masivas está en el orden de Gb de información que sobrepasa la capacidad de memoria de un ordenador de mesa, por lo que podemos usar la función ShortRead::FastqSampler() para obtener una muestra representativa de dicho grupo de datos. Usando la lista de archivos en nuestro directorio de ejemplo fastaqFiles, agregando el raiz del directorio
fastaqFiles[1]## [1] "brca.example.illumina.0.1.fastq"pwd## [1] "C:/Bioinformatica_2020_2_BI0035/data/clase_4/brca/"archivo <- paste0(pwd,fastaqFiles[1])
archivo## [1] "C:/Bioinformatica_2020_2_BI0035/data/clase_4/brca/brca.example.illumina.0.1.fastq"Sacaremos una muestra de 1000 secuencias del primer archivo con la función ShortRead::FastqSampler()
muestra <- ShortRead::FastqSampler(archivo,n = 10000)Sin embargo, lo que hace ShortRead::FastqSampler() es establecer una conexción con cada una de las secuencias de muestras, mas no genera un archivo de secuencias, como lo podemos evidenciar si es que tratamos de imprimir la secuencias en pantalla con la fucnión ShortRead::sread()
# ShortRead::sread(muestra)
# Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘sread’ for signature ‘"FastqSampler"’Nos aparece un rerro indicando que no hay un método “sread” definido para el objeto “FastqSampler”, clase a la que pertence el objeto “muestra”. Para crear un objeto de clase “ShortReadQ” que contenga las secuencias podemos usar la función ShortRead::yield()
m1 <- ShortRead::yield(muestra)
ShortRead::sread(m1)## DNAStringSet object of length 10000:
## width seq
## [1] 101 AGAATGGCATGAACCCGGGAGGCAGAGTTTGC...AACAGAGCAAGACTCCATCTCAAAAAAAGAA
## [2] 101 GGCATGGTAGGGAATATGGGGGAGTGGGAAAG...TGTGAGAATTTGTTGTTGTCGTTGTTGTTGT
## [3] 101 CAGGCAGATGAGTGACAGCAAGAAAACCTGGC...CAATTGGATTTTTATGAAATTTCCAAGTTCC
## [4] 101 TTGCAGTGAGCCAAGATCATACCACGGCACTC...AAAAAAAAAAGGAAAATGAAACTAGAAGAGA
## [5] 101 CTCAGCTCACTGCAACCTCCACCTCCCGGGTT...ATTAAAGGTGCCTGCCACCATGCCAGGCTAA
## ... ... ...
## [9996] 101 ATGTTTTCATCACTGGAACCTATTTCATTAAT...TAATATTGCTTGAGCTGGCTTCTTTAAAAAC
## [9997] 101 TCCTTTGGCCTCTCTGGGTTAGGAAAGTATCC...TTCAAATATAAAACCAGCTATTGGGGGTATG
## [9998] 101 TATGCAAAGTTAGCTTATTTTTCCTGGATGTA...TCTGTTGTCAGGCTGGAGTGCATTGGCATGA
## [9999] 101 AGAAACCTTTGTGAAGCTAGGAATACATTTTC...AGCAGGAGAAAGAAGATGGCTTGTTTTTAGT
## [10000] 101 TCCATGGTGTCAAGTTTCTCTTCAGGAGGAAA...CCATACTGTTTAGCAGGAAACCAGTCTCAGTBiostrings::detail(m1)## class: ShortReadQ
##
## sread:
## DNAStringSet object of length 10000:
## width seq
## [1] 101 AGAATGGCATGAACCCGGGAGGCAGAGTTTGC...AACAGAGCAAGACTCCATCTCAAAAAAAGAA
## [2] 101 GGCATGGTAGGGAATATGGGGGAGTGGGAAAG...TGTGAGAATTTGTTGTTGTCGTTGTTGTTGT
## [3] 101 CAGGCAGATGAGTGACAGCAAGAAAACCTGGC...CAATTGGATTTTTATGAAATTTCCAAGTTCC
## [4] 101 TTGCAGTGAGCCAAGATCATACCACGGCACTC...AAAAAAAAAAGGAAAATGAAACTAGAAGAGA
## [5] 101 CTCAGCTCACTGCAACCTCCACCTCCCGGGTT...ATTAAAGGTGCCTGCCACCATGCCAGGCTAA
## ... ... ...
## [9996] 101 ATGTTTTCATCACTGGAACCTATTTCATTAAT...TAATATTGCTTGAGCTGGCTTCTTTAAAAAC
## [9997] 101 TCCTTTGGCCTCTCTGGGTTAGGAAAGTATCC...TTCAAATATAAAACCAGCTATTGGGGGTATG
## [9998] 101 TATGCAAAGTTAGCTTATTTTTCCTGGATGTA...TCTGTTGTCAGGCTGGAGTGCATTGGCATGA
## [9999] 101 AGAAACCTTTGTGAAGCTAGGAATACATTTTC...AGCAGGAGAAAGAAGATGGCTTGTTTTTAGT
## [10000] 101 TCCATGGTGTCAAGTTTCTCTTCAGGAGGAAA...CCATACTGTTTAGCAGGAAACCAGTCTCAGT
##
## id:
## BStringSet object of length 10000:
## width seq
## [1] 57 Frag_10786 chr17 (Strand + Offset 79673--80234) 562M 101M
## [2] 57 Frag_14504 chr17 (Strand - Offset 15837--16299) 463M 101M
## [3] 55 Frag_15622 chr17 (Strand - Offset 6172--6805) 634M 101M
## [4] 51 Frag_4 chr17 (Strand - Offset 1052--1590) 539M 101M
## [5] 51 Frag_5 chr17 (Strand + Offset 3003--3748) 746M 101M
## ... ... ...
## [9996] 57 Frag_15485 chr17 (Strand + Offset 49034--49541) 508M 101M
## [9997] 57 Frag_10408 chr17 (Strand - Offset 90817--91602) 786M 101M
## [9998] 57 Frag_12999 chr17 (Strand + Offset 70035--70913) 879M 101M
## [9999] 56 Frag_9999 chr17 (Strand - Offset 40946--41411) 466M 101M
## [10000] 57 Frag_13970 chr17 (Strand - Offset 55933--56668) 736M 101M
## class: FastqQuality
## quality:
## BStringSet object of length 10000:
## width seq
## [1] 101 BGDG)GGGGAGGGGGGG?GEF?7FFGFFFGGG...GEF@GDFE??G:F<BGBE=:AECCBCD7@@:
## [2] 101 FGGGFFFGFGDGGCGGG:?FGE:GGGG=FGG;...FC;F7GGCFEGGEEAG?19:AFGGG?@+FDF
## [3] 101 GFCGGGGG9GGAGFGFAE4;BGGFGGBFGGGE...DGG@GGF:GCGGAG;F?CD?GGEEDEDEC7*
## [4] 101 GGEGG?GGGFGEGGG@*F@FFFEGGF>G?GGF...FGEFGCGGFGGGFCCFGGFFF6?E=;6B@FC
## [5] 101 GGFGFEFCFEFG<GGG>82F46GGGEEGC8GG...G0FC(F%FGGC3EFFDAD'9=;16BBGE@<3
## ... ... ...
## [9996] 101 DGGF"GFDCGGGGGEGGGGFG=EGGC=?EGGF...GGG4B3G7GGG=GG@GGGGFF?G=@G?CCEB
## [9997] 101 GGGEDGGFGGGGFGGG@GGGF@GGGFGGG>FF...F5F2GG@CGGGE.EGGGD?BGGFD5A@DE@F
## [9998] 101 GFGCFGGG3BGCAAGGBGFGGGGGGB@BG(FG...'GG4?GAGFGDEA6G@GC:ACF/CEGGEGG6
## [9999] 101 GAGCF1GGGEGGFGEGGBEGGCGEFGFCGGFF...EEGFG8CGEFBGG,GGC@EFFFFG=FGF8FC
## [10000] 101 GGE@FGDF2G?EDGEGFGFGGGEEGCGGFGFG...EBEGGEAGF?GG@FG=AEGFGG66GGEE5=F5. ARCHIVOS FORMATO SAM Y BAM
Este tipo de archivos está dedicado a almacenar la información del alineamiento de las secuencias, la información de la misma secuencia y de su calidad de secuenciación. El formato SAM divide la información en filas y columnas separadas por tabulación, sin embargolas divisiónes, los encabezados y el alineamiento de las secciones lo hace poco amigable para ser interpretado. Una referencia completa de la estructura de estos archivos lo podemos buscar en la página web de Samtools (http://www.htslib.org/doc/sam.html) El formato BAM es una forma comprimida de SAM, que nos ayuda a ahorrar espacio, lo cual es de suma utilidad en los volúmenes de datos que manejamos en secuenciación masiva.
6. FLUJO DE ANALISIS DE DATOS
El flujo de análisis de datos de secuenciación masiva usualmente sigue un patrón uniforme. Lo primero que tenemos que hacer es una revisión de la calidad de la secuenciación para eliminar los lecturas con baja calidad o regiones contaminadas. Luego debemos alinear las secuencias usando programas externos. Una vez alineados, debemos descartar aquellos artefactos de alineamiento. Los pasos finales del análisis incluyen la identificación de variantes, comparaciones estadisticas, anotaciones genómicas y otros análisis específicos requeridos.
6.1 Consideraciones generales del procesamiento de datos de NGS
En general, la mayor dificultad de estos análisis es la memoria. Los datos de secuenciación masiva ocupan un espacio del orden de Gb, por lo que las memorias de las computadoras de mesa o laptops quedan grandemente sobrepasadas. Por ello es que la mayoría de análisis se hace en servidores o clusters. Aún con esta capacidad de computo de 32 o 64 núcleos y memorias RAM de cientos de Gb, estos cálculos pueden durar horas o días. Existen algunas estrategias o recomendaciones para hacer más eficiente el código en R para estos procesos:
- Para tener un ejecución en el código más rápida prefiera usar las funciones de la familia: apply(), mappy(), lapply() que puenden ser implementadas en arreglos, vectores, o listas, en vez de usar loops. Tambien es recomendable que use funciones como Map(), Filter() o Reduce().
- Otra recomendación es que cree estructuras de datos vacias ya sean matrices (matrix()), vectores (vector()) o data frames (data.frame()) para reservar el espacio en la computadora necesario para el procesamiento de los datos.
- Restrinja los procesos a grupos de datos objetivos y haga procesos sobre todo el grupo de datos.
- Puede usar subgrupo de datos para algunos análisis estadísticos, recuerde que R es un lenguje de programación orientado a la estadistica, así que no habrá problemas en implementar o ejecutar pruebas estadisticas inferenciales sobre muestras de la población de datos. Para generar muestaras de datos puede utilizar las funsiones “ShortRead::FastqSampler()” o “Rsamtools::BamSampler()”
- Si la memoria es insuficiente, puede procesar los datos en subgrupos pequeños, usando un loop.
- Para el procesamiento de datos masivos, puede usar clusters, para ello puede reemplazar las funciones apply() por funciones “paralelizadas” como clusterApply(), parApply() o sfApply() de paqueterías como “snow” o “snowfall”
6.2 Verificación de la calidad y escrutinio de las secuencias
El paquete “Biostring” nos ofrece algunas funciones útiles para una explaración rápida de las secuencias. Para genrear tablas de frecuencia absoluta por cada nucleotido podemos usar la función Biostrings::alphabetFrequency sobre las secuencias obtenidas por la función ShortRead::sread()
Biostrings::alphabetFrequency(x = ShortRead::sread(fq),
baseOnly = TRUE,
collapse = TRUE)## A C G T other
## 2620865 2069778 2068976 2623289 29Tambien podemos observar la frecuencia absoluta de nucleótidos de todas las secuencias en cada ciclo con la función ShortRead::alphabetByCycle. Para verlo de una forma mas amigable, podemos graficar estos datos usando la función matplot()
head(ShortRead::alphabetByCycle(ShortRead::sread(fq)))## cycle
## alphabet [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## A 26073 25913 26261 25869 26267 26276 25896 25977 26172 25778 25991
## C 20442 20767 20320 20544 20454 20271 20802 20962 20703 20612 20587
## G 20435 20456 20555 20754 20538 20724 20524 20502 20101 20746 20633
## T 26222 26036 26036 26005 25913 25901 25950 25731 26196 26036 25961
## M 0 0 0 0 0 0 0 0 0 0 0
## R 0 0 0 0 0 0 0 0 0 0 0
## cycle
## alphabet [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## A 26165 26070 25740 26140 25571 26155 26323 25907 26405 26042 26335
## C 20560 20254 20806 20417 20694 20397 20430 20719 20253 20279 20603
## G 20673 20787 20540 20441 20584 20573 20430 20330 20394 20878 20539
## T 25773 26060 26086 26174 26323 26047 25989 26216 26120 25973 25695
## M 0 0 0 0 0 0 0 0 0 0 0
## R 0 0 0 0 0 0 0 0 0 0 0
## cycle
## alphabet [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## A 26025 26178 26533 26159 26472 26277 26159 26106 26061 26015 25940
## C 20633 20596 20526 20438 20638 20352 20430 20528 20628 20641 20887
## G 20433 20487 20350 20799 20289 20534 20572 20631 20567 20620 20451
## T 26081 25911 25763 25776 25772 26008 26009 25906 25915 25895 25893
## M 0 0 0 0 0 0 0 0 0 0 0
## R 0 0 0 0 0 0 0 0 0 0 0
## cycle
## alphabet [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## A 26423 26158 26179 26503 25990 26232 26130 26139 26118 26245 26265
## C 20418 20525 20585 20037 20308 20815 20746 20510 20533 20296 20407
## G 20298 20552 20364 20524 20818 20448 20404 20547 20606 20252 20438
## T 26032 25936 26044 26108 26056 25677 25892 25976 25915 26379 26062
## M 0 0 0 0 0 0 0 0 0 0 0
## R 0 0 0 0 0 0 0 0 0 0 0
## cycle
## alphabet [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55]
## A 26189 26343 25738 26117 25963 26044 26037 25917 25948 26183 26012
## C 20542 20459 20542 20696 20532 20458 20381 20332 20665 20619 20291
## G 20583 20196 20757 20675 20435 20693 20540 20817 20576 20349 20731
## T 25858 26174 26135 25684 26242 25977 26214 26106 25983 26021 26138
## M 0 0 0 0 0 0 0 0 0 0 0
## R 0 0 0 0 0 0 0 0 0 0 0
## cycle
## alphabet [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66]
## A 26023 25919 26107 26194 25838 26060 25897 25976 25776 25976 26128
## C 20507 20823 20549 20588 20533 20957 20403 20663 20689 20666 20434
## G 20551 20232 20267 20316 20727 20467 20449 20566 20699 20576 20511
## T 26091 26198 26249 26074 26074 25688 26423 25967 26008 25954 26098
## M 0 0 0 0 0 0 0 0 0 0 0
## R 0 0 0 0 0 0 0 0 0 0 0
## cycle
## alphabet [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## A 26188 26005 26045 25596 26094 25894 25809 26052 26072 25998 26039
## C 20410 20709 20501 20685 20550 20569 20362 20699 20421 20271 20620
## G 20494 20391 20556 20462 20239 20515 20757 20497 20404 20752 20264
## T 26079 26066 26069 26428 26288 26193 26243 25923 26274 26150 26249
## M 0 0 0 0 0 0 0 0 0 0 0
## R 0 0 0 0 0 0 0 0 0 0 0
## cycle
## alphabet [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
## A 25909 25889 25828 25991 25928 26018 25331 25810 25959 25690 25798
## C 20755 20352 20809 20654 20655 20337 20598 20628 20572 20391 20363
## G 20404 20654 20722 20540 20286 20428 20737 20375 20433 20402 20514
## T 26104 26277 25813 25987 26235 26250 26267 25987 25696 26002 25647
## M 0 0 0 0 0 0 0 0 0 0 0
## R 0 0 0 0 0 0 0 0 0 0 0
## cycle
## alphabet [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99]
## A 25411 25395 25395 25581 25354 25264 25057 25234 25354 25259 25431
## C 20312 20330 20347 20247 19974 20089 20385 20242 19879 20080 20017
## G 20557 20521 20250 20245 20581 20117 20211 20125 20356 20124 20245
## T 25895 25773 25859 25615 25639 25929 25579 25434 25154 25280 25050
## M 0 0 0 0 0 0 0 0 0 0 0
## R 0 0 0 0 0 0 0 0 0 0 0
## cycle
## alphabet [,100] [,101]
## A 25047 25122
## C 20192 20091
## G 20052 19932
## T 25452 25598
## M 0 0
## R 0 0matplot(t(ShortRead::alphabetByCycle(ShortRead::sread(fq))[c("A","C","G","T"),]),
type="l")
En el eje X se define el número del ciclo, y en el eje Y la frecuencia absolta de los nucleótidos.
Otro análisis interesante que podemos hacer es determinar la frecuencia absoluta de los di y tri nucleotidos por cada “read”, usando las funciones Biostrings::dinucleotideFrequency() y Biostrings::trinucleotideFrequency(). Las columnas resultantes de las matrices de estas funciones corresponden a las diferentes combinaciones de di y trinucleótidos, mientas que las filas corresponden a cada uno de los “reads”
# Dinucleotide
head(Biostrings::dinucleotideFrequency(ShortRead::sread(fq)))## AA AC AG AT CA CC CG CT GA GC GG GT TA TC TG TT
## [1,] 6 5 8 5 11 12 1 8 4 6 7 5 2 9 6 5
## [2,] 3 6 6 8 11 15 1 8 7 7 3 2 1 7 9 6
## [3,] 5 4 6 6 6 8 3 11 2 9 2 6 8 8 8 8
## [4,] 23 7 10 3 9 4 1 5 10 5 5 4 2 3 8 1
## [5,] 7 4 8 3 12 14 1 11 1 10 6 3 3 9 5 3
## [6,] 5 4 9 5 6 11 0 6 6 4 6 8 6 4 8 12# Trinucleotide
head(Biostrings::trinucleotideFrequency(ShortRead::sread(fq)))## AAA AAC AAG AAT ACA ACC ACG ACT AGA AGC AGG AGT ATA ATC ATG ATT CAA CAC
## [1,] 2 1 3 0 0 4 0 1 0 2 4 2 1 1 2 1 4 3
## [2,] 1 0 2 0 3 2 0 1 1 3 0 2 0 2 3 3 2 2
## [3,] 1 0 2 2 1 1 1 1 1 1 1 3 1 2 0 3 0 1
## [4,] 17 2 3 1 2 1 1 3 5 2 1 2 1 1 1 0 3 2
## [5,] 1 1 3 1 0 3 0 1 0 4 3 1 0 0 1 2 4 3
## [6,] 1 1 3 0 0 1 0 3 1 2 4 2 2 1 0 2 1 0
## CAG CAT CCA CCC CCG CCT CGA CGC CGG CGT CTA CTC CTG CTT GAA GAC GAG GAT
## [1,] 2 2 5 2 1 4 0 0 1 0 0 5 2 1 0 1 1 2
## [2,] 3 4 7 3 0 5 0 0 1 0 0 4 3 1 0 2 1 4
## [3,] 3 2 2 1 0 4 1 1 1 0 4 3 2 2 2 0 0 0
## [4,] 3 1 3 0 0 1 0 0 1 0 1 2 2 0 3 2 3 1
## [5,] 4 1 5 2 1 6 0 0 1 0 1 7 3 0 0 0 0 1
## [6,] 3 2 3 7 0 1 0 0 0 0 2 0 2 2 2 3 0 1
## GCA GCC GCG GCT GGA GGC GGG GGT GTA GTC GTG GTT TAA TAC TAG TAT TCA TCC
## [1,] 3 2 0 1 1 3 1 2 0 2 2 1 0 0 1 1 3 4
## [2,] 0 6 0 1 1 1 1 0 0 0 2 0 0 1 0 0 1 4
## [3,] 1 4 1 3 0 1 0 1 2 2 1 1 2 3 1 2 2 2
## [4,] 2 2 0 1 1 2 1 1 0 0 4 0 0 1 1 0 2 1
## [5,] 3 5 0 2 1 1 2 2 1 0 1 1 2 0 1 0 4 4
## [6,] 2 1 0 1 3 1 1 1 0 1 4 3 1 0 3 2 1 2
## TCG TCT TGA TGC TGG TGT TTA TTC TTG TTT
## [1,] 0 2 3 1 1 1 1 1 0 2
## [2,] 1 1 5 3 1 0 1 1 1 2
## [3,] 1 3 0 6 0 2 1 1 5 1
## [4,] 0 0 4 1 2 1 0 0 1 0
## [5,] 0 1 0 5 0 0 1 2 0 0
## [6,] 0 1 2 1 1 4 2 2 2 5Podemos sumar la frecuencia absoluta de estas matrices con la función nativa colSums()
# Dinucleotide
colSums(Biostrings::dinucleotideFrequency(ShortRead::sread(fq)))## AA AC AG AT CA CC CG CT GA GC GG
## 842084 478693 683105 591336 667085 551415 146899 683867 544877 474045 551675
## GT TA TC TG TT
## 477751 540746 545182 666859 844112# Trinucleotide
colSums(Biostrings::trinucleotideFrequency(ShortRead::sread(fq)))## AAA AAC AAG AAT ACA ACC ACG ACT AGA AGC AGG
## 350465 126617 169091 187705 166233 123859 30930 152841 191815 146668 186496
## AGT ATA ATC ATG ATT CAA CAC CAG CAT CCA CCC
## 151402 143055 117638 136831 187820 163032 153633 207388 136387 188128 130562
## CCG CCT CGA CGC CGG CGT CTA CTC CTG CTT GAA
## 42297 184996 32503 39874 42174 30694 121172 178053 207214 170460 153409
## GAC GAG GAT GCA GCC GCG GCT GGA GGC GGG GGT
## 89493 179075 117560 137690 144389 40087 147144 146880 145125 129179 124911
## GTA GTC GTG GTT TAA TAC TAG TAT TCA TCC TCG
## 104135 89330 153426 125955 166682 104138 120678 143792 168459 146963 32183
## TCT TGA TGC TGG TGT TTA TTC TTG TTT
## 192062 168223 137675 188363 165928 166996 154548 162670 351373Si nosotros deseamos hacer conteo de grupo de nucleótidos de más de tres grupos, podemos usar la función Biostrings::oligonucleotideFrequency() y podemos definir el número de nucleotidos.
x <- colSums(Biostrings::oligonucleotideFrequency(ShortRead::sread(fq), 5))
matplot(x, type = "l")
a) Revisión de calidad de secuencias simples
Anteriormente ya habíamos creado una muestra de 10000 secuencias para poder trabajar en nuestro ordenador. A ese objeto “ShortReadQ” lo llamamos “m1”
m1## class: ShortReadQ
## length: 10000 reads; width: 101 cyclesBiostrings::detail(m1)## class: ShortReadQ
##
## sread:
## DNAStringSet object of length 10000:
## width seq
## [1] 101 AGAATGGCATGAACCCGGGAGGCAGAGTTTGC...AACAGAGCAAGACTCCATCTCAAAAAAAGAA
## [2] 101 GGCATGGTAGGGAATATGGGGGAGTGGGAAAG...TGTGAGAATTTGTTGTTGTCGTTGTTGTTGT
## [3] 101 CAGGCAGATGAGTGACAGCAAGAAAACCTGGC...CAATTGGATTTTTATGAAATTTCCAAGTTCC
## [4] 101 TTGCAGTGAGCCAAGATCATACCACGGCACTC...AAAAAAAAAAGGAAAATGAAACTAGAAGAGA
## [5] 101 CTCAGCTCACTGCAACCTCCACCTCCCGGGTT...ATTAAAGGTGCCTGCCACCATGCCAGGCTAA
## ... ... ...
## [9996] 101 ATGTTTTCATCACTGGAACCTATTTCATTAAT...TAATATTGCTTGAGCTGGCTTCTTTAAAAAC
## [9997] 101 TCCTTTGGCCTCTCTGGGTTAGGAAAGTATCC...TTCAAATATAAAACCAGCTATTGGGGGTATG
## [9998] 101 TATGCAAAGTTAGCTTATTTTTCCTGGATGTA...TCTGTTGTCAGGCTGGAGTGCATTGGCATGA
## [9999] 101 AGAAACCTTTGTGAAGCTAGGAATACATTTTC...AGCAGGAGAAAGAAGATGGCTTGTTTTTAGT
## [10000] 101 TCCATGGTGTCAAGTTTCTCTTCAGGAGGAAA...CCATACTGTTTAGCAGGAAACCAGTCTCAGT
##
## id:
## BStringSet object of length 10000:
## width seq
## [1] 57 Frag_10786 chr17 (Strand + Offset 79673--80234) 562M 101M
## [2] 57 Frag_14504 chr17 (Strand - Offset 15837--16299) 463M 101M
## [3] 55 Frag_15622 chr17 (Strand - Offset 6172--6805) 634M 101M
## [4] 51 Frag_4 chr17 (Strand - Offset 1052--1590) 539M 101M
## [5] 51 Frag_5 chr17 (Strand + Offset 3003--3748) 746M 101M
## ... ... ...
## [9996] 57 Frag_15485 chr17 (Strand + Offset 49034--49541) 508M 101M
## [9997] 57 Frag_10408 chr17 (Strand - Offset 90817--91602) 786M 101M
## [9998] 57 Frag_12999 chr17 (Strand + Offset 70035--70913) 879M 101M
## [9999] 56 Frag_9999 chr17 (Strand - Offset 40946--41411) 466M 101M
## [10000] 57 Frag_13970 chr17 (Strand - Offset 55933--56668) 736M 101M
## class: FastqQuality
## quality:
## BStringSet object of length 10000:
## width seq
## [1] 101 BGDG)GGGGAGGGGGGG?GEF?7FFGFFFGGG...GEF@GDFE??G:F<BGBE=:AECCBCD7@@:
## [2] 101 FGGGFFFGFGDGGCGGG:?FGE:GGGG=FGG;...FC;F7GGCFEGGEEAG?19:AFGGG?@+FDF
## [3] 101 GFCGGGGG9GGAGFGFAE4;BGGFGGBFGGGE...DGG@GGF:GCGGAG;F?CD?GGEEDEDEC7*
## [4] 101 GGEGG?GGGFGEGGG@*F@FFFEGGF>G?GGF...FGEFGCGGFGGGFCCFGGFFF6?E=;6B@FC
## [5] 101 GGFGFEFCFEFG<GGG>82F46GGGEEGC8GG...G0FC(F%FGGC3EFFDAD'9=;16BBGE@<3
## ... ... ...
## [9996] 101 DGGF"GFDCGGGGGEGGGGFG=EGGC=?EGGF...GGG4B3G7GGG=GG@GGGGFF?G=@G?CCEB
## [9997] 101 GGGEDGGFGGGGFGGG@GGGF@GGGFGGG>FF...F5F2GG@CGGGE.EGGGD?BGGFD5A@DE@F
## [9998] 101 GFGCFGGG3BGCAAGGBGFGGGGGGB@BG(FG...'GG4?GAGFGDEA6G@GC:ACF/CEGGEGG6
## [9999] 101 GAGCF1GGGEGGFGEGGBEGGCGEFGFCGGFF...EEGFG8CGEFBGG,GGC@EFFFFG=FGF8FC
## [10000] 101 GGE@FGDF2G?EDGEGFGFGGGEEGCGGFGFG...EBEGGEAGF?GG@FG=AEGFGG66GGEE5=FPara hacer un análisis de la calidad de las secuencias podemos usar la función ShortRead::qa(). Esta función tiene un pequeño detalle, no se puede generar directamente sobre el objeto “ShortReadQ” generado por la función “ShortRead::yield()”, sino que tenemos que llamar necesarimante a dicha función.
qa1 <- ShortRead::qa(ShortRead::yield(muestra),basename(fastaqFiles[1]))
qa1## class: ShortReadQQA(10)
## QA elements (access with qa[["elt"]]):
## readCounts: data.frame(1 3)
## baseCalls: data.frame(1 5)
## readQualityScore: data.frame(512 4)
## baseQuality: data.frame(95 3)
## alignQuality: data.frame(1 3)
## frequentSequences: data.frame(50 4)
## sequenceDistribution: data.frame(3 4)
## perCycle: list(2)
## baseCall: data.frame(404 4)
## quality: data.frame(3838 5)
## perTile: list(2)
## readCounts: data.frame(0 4)
## medianReadQualityScore: data.frame(0 4)
## adapterContamination: data.frame(1 1)El resultado es un objeto “ShortReadQQA” que es un conjunto de data.frames. Para llamar a uno de ellos solo necesitamos utilizar los dobles corchetes
qa1[["readCounts"]]## read filter aligned
## brca.example.illumina.0.1.fastq 10000 NA NAqa1[["baseCalls"]]## A C G T N
## brca.example.illumina.0.1.fastq 281877 222387 222885 282851 0head(qa1[["readQualityScore"]])## quality density lane type
## 1 0.5257265 4.016532e-06 brca.example.illumina.0.1.fastq read
## 2 0.5979941 1.268440e-05 brca.example.illumina.0.1.fastq read
## 3 0.6702617 3.309138e-05 brca.example.illumina.0.1.fastq read
## 4 0.7425293 7.151045e-05 brca.example.illumina.0.1.fastq read
## 5 0.8147969 1.282166e-04 brca.example.illumina.0.1.fastq read
## 6 0.8870645 1.907635e-04 brca.example.illumina.0.1.fastq readqa1[["baseQuality"]]## score count lane
## 1 0 brca.example.illumina.0.1.fastq
## 2 ! 0 brca.example.illumina.0.1.fastq
## 3 " 322 brca.example.illumina.0.1.fastq
## 4 # 893 brca.example.illumina.0.1.fastq
## 5 $ 2089 brca.example.illumina.0.1.fastq
## 6 % 2191 brca.example.illumina.0.1.fastq
## 7 & 2031 brca.example.illumina.0.1.fastq
## 8 ' 2020 brca.example.illumina.0.1.fastq
## 9 ( 1877 brca.example.illumina.0.1.fastq
## 10 ) 2003 brca.example.illumina.0.1.fastq
## 11 * 2238 brca.example.illumina.0.1.fastq
## 12 + 2321 brca.example.illumina.0.1.fastq
## 13 , 2425 brca.example.illumina.0.1.fastq
## 14 - 2651 brca.example.illumina.0.1.fastq
## 15 . 2997 brca.example.illumina.0.1.fastq
## 16 / 3179 brca.example.illumina.0.1.fastq
## 17 0 3505 brca.example.illumina.0.1.fastq
## 18 1 3794 brca.example.illumina.0.1.fastq
## 19 2 4240 brca.example.illumina.0.1.fastq
## 20 3 4657 brca.example.illumina.0.1.fastq
## 21 4 5069 brca.example.illumina.0.1.fastq
## 22 5 5730 brca.example.illumina.0.1.fastq
## 23 6 6191 brca.example.illumina.0.1.fastq
## 24 7 6856 brca.example.illumina.0.1.fastq
## 25 8 7595 brca.example.illumina.0.1.fastq
## 26 9 8579 brca.example.illumina.0.1.fastq
## 27 : 9447 brca.example.illumina.0.1.fastq
## 28 ; 10590 brca.example.illumina.0.1.fastq
## 29 < 11945 brca.example.illumina.0.1.fastq
## 30 = 13586 brca.example.illumina.0.1.fastq
## 31 > 15071 brca.example.illumina.0.1.fastq
## 32 ? 17536 brca.example.illumina.0.1.fastq
## 33 @ 20236 brca.example.illumina.0.1.fastq
## 34 A 24057 brca.example.illumina.0.1.fastq
## 35 B 29436 brca.example.illumina.0.1.fastq
## 36 C 37340 brca.example.illumina.0.1.fastq
## 37 D 51646 brca.example.illumina.0.1.fastq
## 38 E 81725 brca.example.illumina.0.1.fastq
## 39 F 190207 brca.example.illumina.0.1.fastq
## 40 G 411725 brca.example.illumina.0.1.fastq
## 41 H 0 brca.example.illumina.0.1.fastq
## 42 I 0 brca.example.illumina.0.1.fastq
## 43 J 0 brca.example.illumina.0.1.fastq
## 44 K 0 brca.example.illumina.0.1.fastq
## 45 L 0 brca.example.illumina.0.1.fastq
## 46 M 0 brca.example.illumina.0.1.fastq
## 47 N 0 brca.example.illumina.0.1.fastq
## 48 O 0 brca.example.illumina.0.1.fastq
## 49 P 0 brca.example.illumina.0.1.fastq
## 50 Q 0 brca.example.illumina.0.1.fastq
## 51 R 0 brca.example.illumina.0.1.fastq
## 52 S 0 brca.example.illumina.0.1.fastq
## 53 T 0 brca.example.illumina.0.1.fastq
## 54 U 0 brca.example.illumina.0.1.fastq
## 55 V 0 brca.example.illumina.0.1.fastq
## 56 W 0 brca.example.illumina.0.1.fastq
## 57 X 0 brca.example.illumina.0.1.fastq
## 58 Y 0 brca.example.illumina.0.1.fastq
## 59 Z 0 brca.example.illumina.0.1.fastq
## 60 [ 0 brca.example.illumina.0.1.fastq
## 61 \\ 0 brca.example.illumina.0.1.fastq
## 62 ] 0 brca.example.illumina.0.1.fastq
## 63 ^ 0 brca.example.illumina.0.1.fastq
## 64 _ 0 brca.example.illumina.0.1.fastq
## 65 ` 0 brca.example.illumina.0.1.fastq
## 66 a 0 brca.example.illumina.0.1.fastq
## 67 b 0 brca.example.illumina.0.1.fastq
## 68 c 0 brca.example.illumina.0.1.fastq
## 69 d 0 brca.example.illumina.0.1.fastq
## 70 e 0 brca.example.illumina.0.1.fastq
## 71 f 0 brca.example.illumina.0.1.fastq
## 72 g 0 brca.example.illumina.0.1.fastq
## 73 h 0 brca.example.illumina.0.1.fastq
## 74 i 0 brca.example.illumina.0.1.fastq
## 75 j 0 brca.example.illumina.0.1.fastq
## 76 k 0 brca.example.illumina.0.1.fastq
## 77 l 0 brca.example.illumina.0.1.fastq
## 78 m 0 brca.example.illumina.0.1.fastq
## 79 n 0 brca.example.illumina.0.1.fastq
## 80 o 0 brca.example.illumina.0.1.fastq
## 81 p 0 brca.example.illumina.0.1.fastq
## 82 q 0 brca.example.illumina.0.1.fastq
## 83 r 0 brca.example.illumina.0.1.fastq
## 84 s 0 brca.example.illumina.0.1.fastq
## 85 t 0 brca.example.illumina.0.1.fastq
## 86 u 0 brca.example.illumina.0.1.fastq
## 87 v 0 brca.example.illumina.0.1.fastq
## 88 w 0 brca.example.illumina.0.1.fastq
## 89 x 0 brca.example.illumina.0.1.fastq
## 90 y 0 brca.example.illumina.0.1.fastq
## 91 z 0 brca.example.illumina.0.1.fastq
## 92 { 0 brca.example.illumina.0.1.fastq
## 93 | 0 brca.example.illumina.0.1.fastq
## 94 } 0 brca.example.illumina.0.1.fastq
## 95 ~ 0 brca.example.illumina.0.1.fastqplot(qa1[["baseQuality"]]$count,t="l")
Este objeto tiene mucha información relevante. Una forma de verlo en forma gráfica es usando la función report y proyectandolo en una página web
rurl1 <- report(qa1,dest="quality_example1")## Warning in dir.create(dest, recursive = TRUE): 'quality_example1' already existsbrowseURL(paste0(getwd(), "/quality_example1/index.html"))b) Filtro de secuencias FASTQ por calidad de secuenciación
Uno de los problemas de la secuenciación masiva son las secuencias de baja calidad, que aumentan el ruido de las lecturas. El paquete “ShortReads()” tiene la función nFilter() nos ayuda a eliminar las secuencias de baja calidad. Para ejemplificar esto, filtraremos las lecturas del ejemplo de brca.example.illuina.0.1.fastq.
fq <- ShortRead::readFastq(paste0(dirname(getwd()),"/data/clase_4/brca/",basename(fastaqFiles[1])))
fq.n <- fq[ShortRead::nFilter()(fq)]
fq## class: ShortReadQ
## length: 31713 reads; width: 101 cyclesfq.n## class: ShortReadQ
## length: 31712 reads; width: 101 cyclesTenemos otras tres funciones que nos sirven para “limpiar” nuestras secuenias + trimTails() + trimTailw() + trimEnds()
trimTails()
Esta función lo que va a hacer es eliminar las secuencias en las que ocurran n nucleótidos de baja calidad, definiendo como límite “a”. Esta función tiene 4 parámetros importantes:
- object= donde se define el objeto ShortReadQ que se va a filtrar
- k= el número de nucleótidos seguidos con una calidad definida
- a= donde se define la calidad mínima de la secuencia según el phred score (https://drive5.com/usearch/manual/quality_score.html)
- successive= define si estos deben ser sucesivos o no.
ShortRead::trimTails(object = fq.n,
k = 5,
a = "A",
successive = TRUE)## class: ShortReadQ
## length: 31523 reads; width: 1..101 cyclesfq.trimmed <- ShortRead::trimTails(object = fq.n,
k = 5,
a = "A",
successive = TRUE)
par(mfrow = c(1,2))
plot(density(width(fq.n)), xlim= c(0,130))
plot(density(width(fq.trimmed)), xlim= c(0,130))
Si bien es cierto, ahora tenemos secuencias con mejor calidad, sin embargo, las secuencias con longitudes pequeños son poco deseables para ensamblar. Para filtrar las secuencias con longitudes mayores de 80 nt usamos la función ShortReads::width
fq.trimmed## class: ShortReadQ
## length: 31523 reads; width: 1..101 cyclesfq.trimmed[width(fq.trimmed)>80]## class: ShortReadQ
## length: 29746 reads; width: 81..101 cyclespar(mfrow = c(1,2))
plot(density(width(fq.trimmed)), xlim = c(0,130))
fq.trimmed <- fq.trimmed[width(fq.trimmed)>80]
plot(density(width(fq.trimmed)), xlim = c(0,130))
Para los siguientes pasos, es necesario guardar las secuencias filtradas con la función writeFastq(). Esta función tiene tres parámetros: + object= se define el objeto ShortReadQ que se va a exportar + file= se define la localización del archico donde se va a guardar + compress= se define si se desea comprimir el archivo
ShortRead::writeFastq(object = fq.trimmed,
file = paste0(dirname(getwd()),"/data/clase_4/brca/",basename(fastaqFiles[1]), "_filter1.fastq"),
compress = F
)6.3 Manejo de archivos de alineamientos y variantes genómicas
library(GenomicAlignments)
projectDir <- paste0(dirname(getwd()),"/data/clase_4/brca/")
alnFiles <- list.files(path=projectDir, pattern='*.bam', full.names=T)
aln <- readGAlignments(alnFiles)
aln## GAlignments object with 62985 alignments and 0 metadata columns:
## seqnames strand cigar qwidth start end width
## <Rle> <Rle> <character> <integer> <integer> <integer> <integer>
## [1] chr17 - 101M 101 38450327 38450427 101
## [2] chr17 + 101M 101 38449889 38449989 101
## [3] chr17 + 101M 101 38531147 38531247 101
## [4] chr17 - 101M 101 38531585 38531685 101
## [5] chr17 + 101M 101 38454286 38454386 101
## ... ... ... ... ... ... ... ...
## [62981] chr17 + 101M 101 38515499 38515599 101
## [62982] chr17 - 101M 101 38516199 38516299 101
## [62983] chr17 - 101M 101 38565629 38565729 101
## [62984] chr17 + 101M 101 38564999 38565099 101
## [62985] chr2 + 101M 101 135442574 135442674 101
## njunc
## <integer>
## [1] 0
## [2] 0
## [3] 0
## [4] 0
## [5] 0
## ... ...
## [62981] 0
## [62982] 0
## [62983] 0
## [62984] 0
## [62985] 0
## -------
## seqinfo: 45 sequences from an unspecified genomeComo este archivo tiene 62985 secuencias es posible que produzca un “crash” de la computadora si tratamos de cargar toda la tabla con todas las columnas. Para analizar las columnas de forma indpendiente podemos usar la función table()
Ver a que cromosoma corresponde la secuencia analizada
table(seqnames(aln))##
## chr1 chr2 chr3 chr4 chr5 chr6
## 8 19 11 14 5 4
## chr7 chr8 chr9 chr10 chr11 chr12
## 9 5 5 2 8 10
## chr13 chr14 chr15 chr16 chr17 chr18
## 6 3 4 5 62843 1
## chr19 chr20 chr21 chr22 chrX chrY
## 4 4 0 4 8 3
## chrM chr1_random chr2_random chr3_random chr4_random chr5_random
## 0 0 0 0 0 0
## chr6_random chr7_random chr8_random chr9_random chr10_random chr11_random
## 0 0 0 0 0 0
## chr13_random chr15_random chr16_random chr17_random chr18_random chr19_random
## 0 0 0 0 0 0
## chr21_random chr22_random chrX_random
## 0 0 0Para ver la dircción de las secuencias
table(strand(aln))##
## + - *
## 31489 31496 0table(width(aln))##
## 89 92 93 94 95 96 97 98 99 100 101 102 103
## 2 2 3 2 2 3 29 88 97 229 62274 167 32
## 104 105 108
## 26 28 1head(table(cigar(aln)))##
## 101M 10M1D91M 10M1I90M 10M2D91M 10M2I89M 10M3I88M
## 62259 1 4 2 1 1xtabs(∼seqnames + strand, as.data.frame(aln))## strand
## seqnames + - *
## chr1 6 2 0
## chr2 7 12 0
## chr3 5 6 0
## chr4 6 8 0
## chr5 2 3 0
## chr6 2 2 0
## chr7 6 3 0
## chr8 2 3 0
## chr9 3 2 0
## chr10 1 1 0
## chr11 3 5 0
## chr12 5 5 0
## chr13 3 3 0
## chr14 1 2 0
## chr15 2 2 0
## chr16 3 2 0
## chr17 31418 31425 0
## chr18 1 0 0
## chr19 3 1 0
## chr20 2 2 0
## chr21 0 0 0
## chr22 2 2 0
## chrX 4 4 0
## chrY 2 1 0
## chrM 0 0 0
## chr1_random 0 0 0
## chr2_random 0 0 0
## chr3_random 0 0 0
## chr4_random 0 0 0
## chr5_random 0 0 0
## chr6_random 0 0 0
## chr7_random 0 0 0
## chr8_random 0 0 0
## chr9_random 0 0 0
## chr10_random 0 0 0
## chr11_random 0 0 0
## chr13_random 0 0 0
## chr15_random 0 0 0
## chr16_random 0 0 0
## chr17_random 0 0 0
## chr18_random 0 0 0
## chr19_random 0 0 0
## chr21_random 0 0 0
## chr22_random 0 0 0
## chrX_random 0 0 0xtabs(∼seqnames + width, as.data.frame(aln))## width
## seqnames 89 92 93 94 95 96 97 98 99 100
## chr1 0 0 0 0 0 1 1 0 0 0
## chr2 0 0 0 0 0 0 1 0 3 1
## chr3 0 0 0 0 0 0 1 0 0 0
## chr4 0 0 1 0 0 0 0 0 0 0
## chr5 0 0 0 0 0 0 0 0 0 0
## chr6 0 0 0 0 0 0 1 0 0 0
## chr7 0 0 0 0 0 0 0 0 1 0
## chr8 0 0 0 0 0 0 1 0 1 1
## chr9 0 0 0 0 0 0 1 0 0 0
## chr10 0 0 0 0 0 0 0 0 0 1
## chr11 0 0 0 0 0 0 0 0 1 1
## chr12 0 0 1 0 0 0 0 0 0 1
## chr13 0 0 0 0 0 0 0 1 0 0
## chr14 0 0 0 0 0 0 0 0 0 0
## chr15 0 0 0 0 0 0 0 0 1 0
## chr16 0 0 0 0 0 0 0 0 0 1
## chr17 2 2 1 2 2 2 23 87 89 220
## chr18 0 0 0 0 0 0 0 0 0 0
## chr19 0 0 0 0 0 0 0 0 0 1
## chr20 0 0 0 0 0 0 0 0 0 1
## chr21 0 0 0 0 0 0 0 0 0 0
## chr22 0 0 0 0 0 0 0 0 1 0
## chrX 0 0 0 0 0 0 0 0 0 1
## chrY 0 0 0 0 0 0 0 0 0 0
## chrM 0 0 0 0 0 0 0 0 0 0
## chr1_random 0 0 0 0 0 0 0 0 0 0
## chr2_random 0 0 0 0 0 0 0 0 0 0
## chr3_random 0 0 0 0 0 0 0 0 0 0
## chr4_random 0 0 0 0 0 0 0 0 0 0
## chr5_random 0 0 0 0 0 0 0 0 0 0
## chr6_random 0 0 0 0 0 0 0 0 0 0
## chr7_random 0 0 0 0 0 0 0 0 0 0
## chr8_random 0 0 0 0 0 0 0 0 0 0
## chr9_random 0 0 0 0 0 0 0 0 0 0
## chr10_random 0 0 0 0 0 0 0 0 0 0
## chr11_random 0 0 0 0 0 0 0 0 0 0
## chr13_random 0 0 0 0 0 0 0 0 0 0
## chr15_random 0 0 0 0 0 0 0 0 0 0
## chr16_random 0 0 0 0 0 0 0 0 0 0
## chr17_random 0 0 0 0 0 0 0 0 0 0
## chr18_random 0 0 0 0 0 0 0 0 0 0
## chr19_random 0 0 0 0 0 0 0 0 0 0
## chr21_random 0 0 0 0 0 0 0 0 0 0
## chr22_random 0 0 0 0 0 0 0 0 0 0
## chrX_random 0 0 0 0 0 0 0 0 0 0
## width
## seqnames 101 102 103 104 105 108
## chr1 5 0 0 1 0 0
## chr2 12 2 0 0 0 0
## chr3 9 1 0 0 0 0
## chr4 10 1 2 0 0 0
## chr5 5 0 0 0 0 0
## chr6 3 0 0 0 0 0
## chr7 8 0 0 0 0 0
## chr8 2 0 0 0 0 0
## chr9 4 0 0 0 0 0
## chr10 0 1 0 0 0 0
## chr11 5 1 0 0 0 0
## chr12 8 0 0 0 0 0
## chr13 5 0 0 0 0 0
## chr14 3 0 0 0 0 0
## chr15 2 0 0 0 1 0
## chr16 4 0 0 0 0 0
## chr17 62173 161 28 25 25 1
## chr18 1 0 0 0 0 0
## chr19 3 0 0 0 0 0
## chr20 2 0 0 0 1 0
## chr21 0 0 0 0 0 0
## chr22 3 0 0 0 0 0
## chrX 5 0 2 0 0 0
## chrY 2 0 0 0 1 0
## chrM 0 0 0 0 0 0
## chr1_random 0 0 0 0 0 0
## chr2_random 0 0 0 0 0 0
## chr3_random 0 0 0 0 0 0
## chr4_random 0 0 0 0 0 0
## chr5_random 0 0 0 0 0 0
## chr6_random 0 0 0 0 0 0
## chr7_random 0 0 0 0 0 0
## chr8_random 0 0 0 0 0 0
## chr9_random 0 0 0 0 0 0
## chr10_random 0 0 0 0 0 0
## chr11_random 0 0 0 0 0 0
## chr13_random 0 0 0 0 0 0
## chr15_random 0 0 0 0 0 0
## chr16_random 0 0 0 0 0 0
## chr17_random 0 0 0 0 0 0
## chr18_random 0 0 0 0 0 0
## chr19_random 0 0 0 0 0 0
## chr21_random 0 0 0 0 0 0
## chr22_random 0 0 0 0 0 0
## chrX_random 0 0 0 0 0 0