How to deal with a file.info and missing value error in R

user4cgI 0
Viewed 77

I am trying to use the following example code on my binary alignment (.bam) file. The example code is available at https://rdrr.io/bioc/ATACseqQC/man/TSSEscore.html

library(GenomicRanges)
bamfile <- system.file("extdata", "GL1.bam", 
                   package="ATACseqQC", mustWork=TRUE)
gal1 <- readBamFile(bamFile=bamfile, tag=character(0), 
                which=GRanges("chr1", IRanges(1, 1e6)), 
                asMates=FALSE)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txs <- transcripts(TxDb.Hsapiens.UCSC.hg19.knownGene)
tsse <- TSSEscore(gal1, txs)

I am getting the following error of missing value TRUE/FALSE when I am reading my bamfile

bamfile < readBamFile("extdata","RRrep1_rmDupChrC_sorted.bam", package="ATACseqQC", mustWork="TRUE")

Seen below:

Error in if (file.size(bamFile) > 1e+08 && !bigFile && interactive()) { : 
  missing value where TRUE/FALSE needed

To fix this error I added the following if condition before reading my bam file:

RRrep1_rmDupChrC_sorted.bam = TRUE
if(file.size(RRrep1_rmDupChrC_sorted.bam) > 1e+08 && !bigFile && interactive()) {RRrep1_rmDupChrC_sorted.bam=FALSE} else {if(RRrep1_rmDupChrC_sorted.bam) {RRrep1_rmDupChrC_sorted.bam}}

At the above step I get this error:

Error in file.info(..., extra_cols = FALSE) : invalid filename argument

Below is a part of my input bamfile using dput.

dput(head(mybamfile[100]))   
new("GAlignments", seqnames = new("Rle", values = 
structure(1L, levels = c("chr_1A",

"chr_1B", "chr_1C", "chr_1D", "chr_2A", "chr_2B", "chr_2C", "chr_2D", "chr_3A", "chr_3B", "chr_3C", "chr_3D", "chr_4A", "chr_4B", "chr_4C", "chr_4D", "chr_5A", "chr_5B", "chr_5C", "chr_5D", "chr_6A", "chr_6B", "chr_6C", "chr_6D", "chr_7A", "chr_7B", "chr_7C", "chr_7D", "contig_01", "contig_02", "contig_03", "contig_04", "contig_05", "contig_06", "contig_07", "contig_08", "contig_09", "contig_10", "contig_11", "contig_12", "contig_13", "contig_14", "contig_15", "contig_16", "contig_17", "contig_18", "contig_19", "contig_20", "contig_21" ), class = "factor"), lengths = 1L, elementMetadata = NULL, metadata = list()), start = 86L, cigar = "89M", strand = new("Rle", values = structure(1L, levels = c("+", "-", "*"), class = "factor"), lengths = 1L, elementMetadata = NULL, metadata = list()), NAMES = NULL, elementMetadata = new("DFrame", rownames = NULL, nrows = 1L, elementType = "ANY", elementMetadata = NULL, metadata = list(), listData = list(qname = "A01488:93:H2YM7DRX2:2:2176:2618:14826", flag = 99L, mapq = 1L, isize = 694L, seq = new("DNAStringSet", pool = new("SharedRaw_Pool", xp_list = list(<pointer: 0x0>), .link_to_cached_object_list = list()), ranges = new("GroupedIRanges", group = 1L, start = 9097L, width = 89L, NAMES = NULL, elementType = "ANY", elementMetadata = NULL, metadata = list()), elementType = "DNAString", elementMetadata = NULL, metadata = list()), qual = new("PhredQuality", pool = new("SharedRaw_Pool", xp_list = list(<pointer: 0x0>), .link_to_cached_object_list = list()), ranges = new("GroupedIRanges", group = 1L, start = 9097L, width = 89L, NAMES = NULL, elementType = "ANY", elementMetadata = NULL, metadata = list()), elementType = "BString", elementMetadata = NULL, metadata = list()), mrnm = structure(1L, levels = c("chr_1A", "chr_1B", "chr_1C", "chr_1D", "chr_2A", "chr_2B", "chr_2C", "chr_2D", "chr_3A", "chr_3B", "chr_3C", "chr_3D", "chr_4A", "chr_4B", "chr_4C", "chr_4D", "chr_5A", "chr_5B", "chr_5C", "chr_5D", "chr_6A", "chr_6B", "chr_6C", "chr_6D", "chr_7A", "chr_7B", "chr_7C", "chr_7D", "contig_01", "contig_02", "contig_03", "contig_04", "contig_05", "contig_06", "contig_07", "contig_08", "contig_09", "contig_10", "contig_11", "contig_12", "contig_13", "contig_14", "contig_15", "contig_16", "contig_17", "contig_18", "contig_19", "contig_20", "contig_21"), class = "factor"))), seqinfo = new("Seqinfo", seqnames = c("chr_1A", "chr_1B", "chr_1C", "chr_1D", "chr_2A", "chr_2B", "chr_2C", "chr_2D", "chr_3A", "chr_3B", "chr_3C", "chr_3D", "chr_4A", "chr_4B", "chr_4C", "chr_4D", "chr_5A", "chr_5B", "chr_5C", "chr_5D", "chr_6A", "chr_6B", "chr_6C", "chr_6D", "chr_7A", "chr_7B", "chr_7C", "chr_7D", "contig_01", "contig_02", "contig_03", "contig_04", "contig_05", "contig_06", "contig_07", "contig_08", "contig_09", "contig_10", "contig_11", "contig_12", "contig_13", "contig_14", "contig_15", "contig_16", "contig_17", "contig_18", "contig_19", "contig_20", "contig_21"), seqlengths = c(22855991L, 25931632L, 25722849L, 26286749L, 27491849L, 27417387L, 27223255L, 28082264L, 31413151L, 31256976L, 31136822L, 30162378L, 26690241L, 30438350L, 26803992L, 24913570L, 26406975L, 29153715L, 27044276L, 26110821L, 34436226L, 35982586L, 35269636L, 32927424L, 22554989L, 23234183L, 22854303L, 22595168L, 603133L, 557894L, 540709L, 344755L, 268859L, 188224L, 173710L, 167197L, 161305L, 158584L, 149526L, 133680L, 133252L, 130241L, 126950L, 75620L, 73032L, 51670L, 49619L, 32731L, 24925L), is_circular = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), genome = c(NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_)), elementType = "ANY", metadata = list(header = c("@HD\tVN:1.6\tSO:coordinate", "@SQ\tSN:chr_1A\tLN:22855991", "@SQ\tSN:chr_1B\tLN:25931632", "@SQ\tSN:chr_1C\tLN:25722849", "@SQ\tSN:chr_1D\tLN:26286749", "@SQ\tSN:chr_2A\tLN:27491849", "@SQ\tSN:chr_2B\tLN:27417387", "@SQ\tSN:chr_2C\tLN:27223255", "@SQ\tSN:chr_2D\tLN:28082264", "@SQ\tSN:chr_3A\tLN:31413151", "@SQ\tSN:chr_3B\tLN:31256976", "@SQ\tSN:chr_3C\tLN:31136822", "@SQ\tSN:chr_3D\tLN:30162378", "@SQ\tSN:chr_4A\tLN:26690241", "@SQ\tSN:chr_4B\tLN:30438350", "@SQ\tSN:chr_4C\tLN:26803992", "@SQ\tSN:chr_4D\tLN:24913570", "@SQ\tSN:chr_5A\tLN:26406975", "@SQ\tSN:chr_5B\tLN:29153715", "@SQ\tSN:chr_5C\tLN:27044276", "@SQ\tSN:chr_5D\tLN:26110821", "@SQ\tSN:chr_6A\tLN:34436226", "@SQ\tSN:chr_6B\tLN:35982586", "@SQ\tSN:chr_6C\tLN:35269636", "@SQ\tSN:chr_6D\tLN:32927424", "@SQ\tSN:chr_7A\tLN:22554989", "@SQ\tSN:chr_7B\tLN:23234183", "@SQ\tSN:chr_7C\tLN:22854303", "@SQ\tSN:chr_7D\tLN:22595168", "@SQ\tSN:contig_01\tLN:603133", "@SQ\tSN:contig_02\tLN:557894", "@SQ\tSN:contig_03\tLN:540709", "@SQ\tSN:contig_04\tLN:344755", "@SQ\tSN:contig_05\tLN:268859", "@SQ\tSN:contig_06\tLN:188224", "@SQ\tSN:contig_07\tLN:173710", "@SQ\tSN:contig_08\tLN:167197", "@SQ\tSN:contig_09\tLN:161305", "@SQ\tSN:contig_10\tLN:158584", "@SQ\tSN:contig_11\tLN:149526", "@SQ\tSN:contig_12\tLN:133680", "@SQ\tSN:contig_13\tLN:133252", "@SQ\tSN:contig_14\tLN:130241", "@SQ\tSN:contig_15\tLN:126950", "@SQ\tSN:contig_16\tLN:75620", "@SQ\tSN:contig_17\tLN:73032", "@SQ\tSN:contig_18\tLN:51670", "@SQ\tSN:contig_19\tLN:49619", "@SQ\tSN:contig_20\tLN:32731", "@SQ\tSN:contig_21\tLN:24925", "@PG\tID:bowtie2\tPN:bowtie2\tVN:2.4.4\tCL:"/mnt/home/agarw118/anaconda3/envs/bowtie2/bin/bowtie2-align-s --wrapper basic-0 -x farr1_chrc -p 8 --very-sensitive -X 1000 -k 10 --no-mixed --dovetail -S RRrep1_align.bam -1 /mnt/gs21/scratch/agarw118/strawberry_almond_atacseq/strawberry/RRrep1_r1.trim.fq.gz -2 /mnt/gs21/scratch/agarw118/strawberry_almond_atacseq/strawberry/RRrep1_r2.trim.fq.gz"", "@PG\tID:samtools\tPN:samtools\tPP:bowtie2\tVN:1.14\tCL:samtools sort -@ 15 -o RRrep1_align_sorted.bam RRrep1_align.bam", "@PG\tID:MarkDuplicates\tVN:2.26.10\tCL:MarkDuplicates INPUT=[RRrep1_align_sorted.bam] OUTPUT=RRrep1_rmDup.bam METRICS_FILE=RRrep1_rmDup.bam_metrics.txt REMOVE_DUPLICATES=true TMP_DIR=[tmp] VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP=50000 MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 SORTING_COLLECTION_SIZE_RATIO=0.25 TAG_DUPLICATE_SET_MEMBERS=false REMOVE_SEQUENCING_DUPLICATES=false TAGGING_POLICY=DontTag CLEAR_DT=true DUPLEX_UMI=false ADD_PG_TAG_TO_READS=true ASSUME_SORTED=false DUPLICATE_SCORING_STRATEGY=SUM_OF_BASE_QUALITIES PROGRAM_RECORD_ID=MarkDuplicates PROGRAM_GROUP_NAME=MarkDuplicates READ_NAME_REGEX=<optimized capture of last three ':' separated fields as numeric values> OPTICAL_DUPLICATE_PIXEL_DISTANCE=100 MAX_OPTICAL_DUPLICATE_SET_SIZE=300000 VERBOSITY=INFO QUIET=false COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json USE_JDK_DEFLATER=false USE_JDK_INFLATER=false\tPN:MarkDuplicates", "@PG\tID:samtools.1\tPN:samtools\tPP:samtools\tVN:1.14\tCL:samtools view -@ 35 -h RRrep1_rmDup.bam", "@PG\tID:samtools.2\tPN:samtools\tPP:samtools.1\tVN:1.14\tCL:samtools sort -@ 35 -o RRrep1_rmDupChrC_sorted.bam RRrep1_rmDupChrC.bam" )))

r
1 Answers

Just as a side note, installing packages for this was an annoying process compared to the normal way in R. You can skip the first two steps as I'm sure you already did them, but the rest of this code worked for me for the first part you ran:

#### Setup BiocManager Package (Skip if Done) ####
if (!require("BiocManager", 
             quietly = TRUE))
  install.packages("BiocManager")

#### Install Necessary Packages (Skip if Done) ####
BiocManager::install(c("ATACseqQC",
                       "TxDb.Hsapiens.UCSC.hg19.knownGene",
                       "GenomicAlignments"))

#### Load Libraries ####
library(GenomicAlignments)
library(GenomicRanges)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(ATACseqQC)

#### Read in System File ####
bamfile <- system.file("extdata", 
                       "GL1.bam", 
                       package="ATACseqQC",
                       mustWork=TRUE)
bamfile

gal1 <- readBamFile(bamFile=bamfile, 
                    tag=character(0), 
                    which=GRanges("chr1", 
                                  IRanges(1, 
                                          1e6)), 
                    asMates=FALSE)

#### Finally Run These ####
txs <- transcripts(TxDb.Hsapiens.UCSC.hg19.knownGene)
tsse <- TSSEscore(gal1, 
                  txs)
tsse

Which gave me this:

$values
         1          2          3          4          5          6 
 2.6972021  3.6455679  2.1226520  3.3067057  5.4685081  4.2493540 
         7          8          9         10         11         12 
 4.6026253  8.5089197  8.2507251 10.9537203  6.0820541  9.6279959 
        13         14         15         16         17         18 
11.4431240  1.0558437  6.7524966  3.0343532  1.0194494  3.9691559 
        19         20 
 0.9155013  1.1102246 

$TSSEscore
[1] 11.44312

For your file, I saved your dput as so.bam. When I ran that object, it gave me this:

GAlignments object with 0 alignments and 0 metadata columns:
   seqnames strand       cigar    qwidth     start       end     width
      <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
       njunc
   <integer>
  -------
  seqinfo: no sequences

So from my end, the bam file seemed to be loaded fine, but perhaps you had a different issue. Judging from this output, it looks as though the file is empty, but I'm not sure if that is a different problem from what you are having.

Shawn Hemelstrand0
Related