Job script

 htseq-count \
    /export/galaxy-central/database/files/000/dataset_41.dat /export/galaxy-central/database/files/000/dataset_42.dat /export/galaxy-central/database/files/000/dataset_43.dat /export/galaxy-central/database/files/000/dataset_44.dat \
    /export/galaxy-central/database/files/000/dataset_45.dat \
    -f bam \
    -r name \
    -s yes \
    -a 10 \
    -t exon \
    -i ID \
    -m union > htseq-counts-raw.txt

  grep -v '__no_feature\|__ambiguous\|__too_low_aQual\|__not_aligned\|__alignment_not_unique' htseq-counts-raw.txt > counts.txt

Counts

Write data into a CSV file.

count_data = read.table(paste0(opt$X_d, '/counts.txt'), row.names = 1)
sample_names = trimws(strsplit(opt$X_B, ',')[[1]])
colnames(count_data) = rep(sample_names, length = ncol(count_data))
# kable(head(count_data))

# modify column names
count_data = data.frame(feature_id = rownames(count_data), count_data)
write.csv(count_data, 
          file = paste0(Sys.getenv('REPORT_FILES_PATH'), '/count_data.csv'),
          quote = FALSE, row.names = FALSE)

Display the top 1000 rows with largest average counts.

# Sort count table by count average
rownames(count_data) = count_data$feature_id
count_data = count_data[, -1]
sorted_ct_table = count_data[order(rowMeans(count_data), decreasing = TRUE), ]
DT::datatable(head(sorted_ct_table, 1000))

Outputs