I know this question is answered many times for different scenarios and I know in the vignette, there is a detailed explanation on that. However, I'm unable to understand the logic behind and genarate a logical/working design matrix for my analysis.

What I want to do is: identify DE genes for KRAS WT and MUT group by controlling the patient_id of tumor-organoid pairs together with RIN values and type of sample (organoid or tumor).

This is my coldata:

# coldata

coldata <- data.frame(
  Patient_ID = c("pat1", "pat1", "pat2", "pat2", "pat3", "pat3", "pat4", "pat4", "pat5", "pat5"),
  Type = c("T", "O", "T", "O", "T", "O", "T", "O", "T", "O"),
  KRAS_Status = c("WT", "WT", "WT", "WT", "MUT", "MUT", "MUT", "MUT", "WT", "WT"),
  RIN = c(5.2, 9.6, 5.6, 9.4, 3.7, 9.6, 2.5, 9.2, 5.2, 8.1),
  Scaled_RIN = c(-0.5311567, 1.4826618, -0.3480823, 1.3911247, -1.2176857, 1.4826618, -1.7669090, 1.2995874, -0.5311567, 0.7961328),
  row.names = c("tum1", "org1", "tum2", "org2", "tum3", "org3", "tum4", "org4", "tum5", "org5")
)

#This is the design that I want:

dds <- DESeqDataSetFromMatrix(countData = round(raw_counts),
                              colData = coldata,
                              design = ~ Patient_ID + KRAS_Status + Type + Scaled_RIN)

I am getting Model matrix not full rank error and will appreciate if you can help with the desig matrix.

You cannot test for KRAS while blocking for patient because patients 1, 2 and 5 are nested with WT and 3-4 are nested with MUT.

What you technically can do in DESeq2 is ~KRAS_Status + Type (I ignore RIN here for simplicity). If you want to account for the correlation between T and O per patient you can use duplicateCorrelation from limma, but then have to go along with limma-voom or limma-trend for the testing.