error fread/samtools

[ltosti]

Hi there,

I performed DrocNc-seq and zUMI seems the ideal tool to analyze my intron-reach data. Unfortunately I have been trying to run zUMI for a couple of days but I get this error and I'm not sure what is causing it; here is the log:

samtools view: writing to standard output failed: Broken pipe
samtools view: error closing standard output: -1
Taking input= as a system command ('samtools view /mnt//lt88/zUMI_AFES1.filtered.tagged.bam | cut -f10 | head -n 100') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.

EXITING: FATAL INPUT ERROR: unrecoginzed parameter name "readFilesType" in input "Command-Line-Initial"
SOLUTION: use correct parameter name (check the manual)

Oct 05 15:28:02 ...... FATAL ERROR, exiting
Error in file(file, "rt") : cannot open the connection
Calls: .runFeatureCount -> <Anonymous> -> read.delim -> read.table -> file
In addition: Warning message:
In file(file, "rt") :
  cannot open file './.Rsubread_featureCounts_pid57353': No such file or directory
Execution halted
Error in gzfile(file, "rb") : cannot open the connection
Calls: readRDS -> gzfile
In addition: Warning message:
In gzfile(file, "rb") :
  cannot open compressed file '/mnt/lt88/zUMI/zUMIs_output/expression/zUMI_AFES1.dgecounts.rds', probable reason 'No such file or directory'
Execution halted

Any idea? Thank you for your help!

[sdparekh]

Hi, to me it looks like you have read and write permission issues. Can you also share your YAML file and the command you are using.

Have you successfully installed all the dependencies as described here: https://github.com/sdparekh/zUMIs/wiki/Installation

[ltosti]

Hi, Thank you for your quick reply. I'm working on a cluster, CentOs 7.5.1804 with TORQUE job scheduler and all the dependencies were installed in a conda environment. This is my YAML file:

project: zUMI_AFES1
sequence_files:
  file1:
    name: /mnt/lt88/fastq/AFES448F1_S1_L001_R1_001.fastq.gz
    base_definition:
    - BC(1-16)
    - UMI(17-26)
  file2:
    name: /mnt/lt88/fastq/AFES448F1_S1_L001_R2_001.fastq.gz
    base_definition: cDNA(49)
reference:
  STAR_index: /mnt/lt88/STAR_transcriptome_zUMI
  GTF_file: /mnt/lt88/refdata-cellranger-hg19-1.2.0/genes/genes.gtf
  additional_STAR_params: ''
  additional_files: ~
out_dir: /mnt/lt88/zUMI
num_threads: 8
mem_limit: 0
filter_cutoffs:
  BC_filter:
    num_bases: 1
    phred: 20
  UMI_filter:
    num_bases: 1
    phred: 20
barcodes:
  barcode_num: ~
  barcode_file: /mnt/lt88/zUMI_bc/barcodes.txt
  BarcodeBinning: 0
  nReadsperCell: 100
counting_opts:
  introns: yes
  downsampling: '0'
  strand: 0
  Ham_Dist: 0
  velocyto: no
  primaryHit: no
  twoPass: yes
make_stats: yes
which_Stage: Filtering

[sdparekh]

Can you send me the output of ls -lhR /mnt/lt88/zUMI/. I am curious if the filtering was finished. Another question is if you gave the paths to samtools and STAR or are they in the PATH variable called as samtools and STAR.

[ltosti]

Hi,

I started another run but I forgot to change the out_dir, so I guess files will be over-written. Sorry, my mistake. Anyway, I am just trying to use another version of samtools; as soon as it finishes, if I get the same error (I assume I will) I'll send you the output of ls -lhR /mnt/lt88/zUMI/. About samtools and STAR they are in the PATH called as samtools and STAR.

[ltosti]

Ok, I got the same error. Here is the output of ls -lhR /mnt/lt88/zUMI/:

zUMI:
total 5.1G
-rwxrwx--- 1 lm15 nobody 142M Oct  5 19:59 zUMI_AFES1.BCstats.txt
-rwxrwx--- 1 lm15 nobody 4.9G Oct  5 19:59 zUMI_AFES1.filtered.tagged.bam
-rwxrwx--- 1 lm15 nobody 1.5K Oct  5 19:59 zUMI_AFES1.postmap.yaml
drwxrws--- 4 lm15 nobody 4.0K Oct  5 20:01 zUMIs_output
-rwxrwx--- 1 lm15 nobody  314 Oct  5 18:26 zUMIs_runlog.txt

zUMI/zUMIs_output:
total 2.3M
drwxrws--- 2 lm15 nobody 4.0K Oct  5 18:25 expression
drwxrws--- 2 lm15 nobody 4.0K Oct  5 18:25 stats
-rwxrwx--- 1 lm15 nobody 2.2M Oct  5 20:01 zUMI_AFES1kept_barcodes.txt

zUMI/zUMIs_output/expression:
total 0

zUMI/zUMIs_output/stats:
total 0

Any idea? Thanks!

[cziegenhain]

Hey,

Sorry that you are having troubles. It seems like the filtering step finishes without issues. However, seems like the issue occurs in the little R piece that checks the data prior to mapping. It requires to call samtools from R and will try to write into the default R temporary directory.

Could you check which temporary directory your R within TORQUE defaults to? tempdir(). Maybe it's a permission issue there?

Best, Christoph

[ltosti]

The tempdir is a directory where I have full access; I agree the problem is is in the R piece, just need to sort out which step!

[cziegenhain]

I have this line in mind: https://github.com/sdparekh/zUMIs/blob/master/zUMIs-mapping.R#L57

Not sure why the samtools view there would fail though... Would you like to manually try in R test <- data.table::fread("samtools view zUMI_AFES1.filtered.tagged.bam | cut -f10 | head -n 100",stringsAsFactors = F,data.table = T, header = F)

[ltosti]

That works fine:

> test <- data.table::fread("samtools view zUMI_AFES1.filtered.tagged.bam | cut -f10 | head -n 100",stringsAsFactors = F,data.table = T, header = F)
> test
     V1
  1:  *
  2:  *
  3:  *
  4:  *
  5:  *
  6:  *
  7:  *
.
.
.
100: *

Is it the output that you would expect?

[cziegenhain]

No that looks actually odd, can we see the first few lines of the bam file?

[ltosti]

samtools view zUMI_AFES1.filtered.tagged.bam | head -2
NB501764:588:HL52TBGX7:1:11101:4432:1065    4   *   0   0   *   *   0   0   *   *   BC:Z:CAACCAATCTCGTATT   UB:Z:TGTTATTTTT
NB501764:588:HL52TBGX7:1:11101:22178:1065   4   *   0   0   *   *   0   0   *   *   BC:Z:GCTTGAAGTCGCTTTC   UB:Z:CTAATTTAAG

[cziegenhain]

Ah! It seems like the cDNA sequence is missing. Sorry I didn't spot that error in the YAML you had posted, it should be cDNA(1-49)

Could you please retry with that corrected? If the error still occurs, would you mind sending us 1mio reads to test on?

[ltosti]

Ah, that would make sense, I will try now. Just for you to know, I used the online shiny app to generate the YAML file and I didn't modify it...I wonder if there is a possible error there? Anyway, I'll keep you posted, thank you for all your help!!

[ltosti]

Hi, I fixed the cDNA range in the YAML file, but unfortunately I just got the same error:

samtools view: writing to standard output failed: Broken pipe
samtools view: error closing standard output: -1
Taking input= as a system command ('samtools view 'samtools view /mnt/lt88/zUMI_AFES1.filtered.tagged.bam | cut -f10 | head -n 100') and a variable has been used in the expression passed to `input=`. Please use fread(cmd=...). There is a security concern if you are creating an app, and the app could have a malicious user, and the app is not running in a secure envionment; e.g. the app is running as root. Please read item 5 in the NEWS file for v1.11.6 for more information and for the option to suppress this message.

EXITING: FATAL INPUT ERROR: unrecoginzed parameter name "readFilesType" in input "Command-Line-Initial"
SOLUTION: use correct parameter name (check the manual)

I tried again to do it in R and this time looks correct?

> test <- data.table::fread("samtools view zUMI_AFES1.filtered.tagged.bam | cut -f10 | head -n 100",stringsAsFactors = F,data.table = T, header = F)
> head(test,2)
                                                  V1
1: GGGGGGGGGAGGTAGAGAGGGGGGGGGATCAACACACACCCAACCCCCC
2: GCAAGAACAGGCTTGTCTAAAGGGAGAGGTACCACCCCCACGACAAGCA

[sdparekh]

Hi, do you mind sharing some of your data? Maybe subset of 100,000 reads? you can get my email address from our paper https://doi.org/10.1093/gigascience/giy059 (First author)

[ltosti]

Hi @sdparekh @cziegenhain , I was using STAR 2.5.3a and as you guessed that was the problem. After updating, the pipeline ran smoothly. Thank you for your time and help!

[cziegenhain]

Glad to hear that! Let us know if we can help further.