hagarelsayed / Ngs_2nd_Abstract

Second Abstract approved by Dr Tamer for RNA_Seq Differential Expression Analysis
0 stars 0 forks source link

Troubleshooting the deseq error #6

Open hagarelsayed opened 4 years ago

hagarelsayed commented 4 years ago

Troubleshooting the deseq error :

An error popped up inferring that;

Error in parametricDispersionFit(means, disps) :
Parametric dispersion fit failed. Try a local fit and/or a pooled estimation.

err

Aligning to our original index (not a subset like in the past example) :

There was a direction towards changing the reference to our original reference that gave high alignment rate with the dataset so that the feature count will result in a count matrix containing different versatile values for deseq to work on.

INDEX=gencode.v33.transcripts

###Other variables are the same and defined previously 
RUNLOG=runlog.txt
READS_DIR=~/workdir/sample_data/renamed/
mkdir bam

for SAMPLE in UNT; do     for REPLICATE in 12 16 20;     do         R1=$READS
_DIR/${SAMPLE}_Rep${REPLICATE}*pass_1.fastq;         R2=$READS_DIR/${SAMPLE}_Rep${REPLICATE}*pass_2.fastq;         BAM=bam/${SAMPLE}_${REPLICATE
}.bam;          hisat2 $INDEX -1 $R1 -2 $R2 | samtools sort > $BAM;         samtools index $BAM;     done; done

Results of the Alignment score for the untreated condition 500 reads; of these: 500 (100.00%) were paired; of these: 71 (14.20%) aligned concordantly 0 times 108 (21.60%) aligned concordantly exactly 1 time 321 (64.20%) aligned concordantly >1 times

71 pairs aligned concordantly 0 times; of these:
  4 (5.63%) aligned discordantly 1 time
----
67 pairs aligned 0 times concordantly or discordantly; of these:
  134 mates make up the pairs; of these:
    112 (83.58%) aligned 0 times
    9 (6.72%) aligned exactly 1 time
    13 (9.70%) aligned >1 times

88.80% overall alignment rate 500 reads; of these: 500 (100.00%) were paired; of these: 75 (15.00%) aligned concordantly 0 times 121 (24.20%) aligned concordantly exactly 1 time 304 (60.80%) aligned concordantly >1 times

75 pairs aligned concordantly 0 times; of these:
  2 (2.67%) aligned discordantly 1 time
----
73 pairs aligned 0 times concordantly or discordantly; of these:
  146 mates make up the pairs; of these:
    132 (90.41%) aligned 0 times
    4 (2.74%) aligned exactly 1 time
    10 (6.85%) aligned >1 times

86.80% overall alignment rate 500 reads; of these: 500 (100.00%) were paired; of these: 86 (17.20%) aligned concordantly 0 times 99 (19.80%) aligned concordantly exactly 1 time 315 (63.00%) aligned concordantly >1 times

86 pairs aligned concordantly 0 times; of these:
  1 (1.16%) aligned discordantly 1 time
----
85 pairs aligned 0 times concordantly or discordantly; of these:
  170 mates make up the pairs; of these:
    147 (86.47%) aligned 0 times
    8 (4.71%) aligned exactly 1 time
    15 (8.82%) aligned >1 times

85.30% overall alignment rate

Same step for the other group

for SAMPLE in UNT; do     for REPLICATE in 12 16 20;     do         R1=$READS
_DIR/${SAMPLE}_Rep${REPLICATE}*pass_1.fastq;         R2=$READS_DIR/${SAMPLE}_Rep${REPLICATE}*pass_2.fastq;         BAM=bam/${SAMPLE}_${REPLICATE
}.bam;          hisat2 $INDEX -1 $R1 -2 $R2 | samtools sort > $BAM;         samtools index $BAM;     done; done

Alignment summary results :+1: 500 reads; of these: 500 (100.00%) were paired; of these: 87 (17.40%) aligned concordantly 0 times 111 (22.20%) aligned concordantly exactly 1 time 302 (60.40%) aligned concordantly >1 times

87 pairs aligned concordantly 0 times; of these:
  0 (0.00%) aligned discordantly 1 time
----
87 pairs aligned 0 times concordantly or discordantly; of these:
  174 mates make up the pairs; of these:
    143 (82.18%) aligned 0 times
    6 (3.45%) aligned exactly 1 time
    25 (14.37%) aligned >1 times

85.70% overall alignment rate

err2

Quantification

GTF=~/workdir/hisat_align/hisatIndex/proj_index/gencode.v33.annotation.gtf
featureCounts -a $GTF -g gene_name -o counts.txt  bam/UNT*.bam  bam/TTT*.bam

Result of Quantification

Differential Expression:

cat simple_counts.txt | Rscript deseq1.r 3x3 > results_deseq1.tsv Another Error popped up ; 33

The table of the count

cat results_deseq1.tsv | awk ' $8 < 0.05 { print $0 }' > filtered_results_deseq1.tsv
cat filtered_results_deseq1.tsv | Rscript draw-heatmap.r > hisat_output.pdf