The dataset used in this course is from a local researcher Weici Zhang generated in the second half of 2020.

I have manipulated the original samples to:

1. combine samples into a single sample.
2. reduce the number of reads for the example dataset to a more managable number for the course.

Data Setup

Let’s set up a project directory for the analysis, and talk a bit about project philosophy.

1. First, create a directory for your user and the example project in the workshop directory:

cd
mkdir -p /share/workshop/intro_scrnaseq/$USER/scrnaseq_example

2a. Next, go into that directory, create a raw data directory (we are going to call this 00-RawData) and cd into that directory. Lets then create symbolic links to the fastq files that contains the raw read data.

cd /share/workshop/intro_scrnaseq/$USER/scrnaseq_example
mkdir 00-RawData
cd 00-RawData/
ln -s /share/workshop/intro_scrnaseq/raw_data/00-SmRawData/PBMC2sm_S18* .

This directory now contains the reads for each “sample” (in this case just 1).

2b. lets create a sample sheet for the project, store sample names in a file called samples.txt

echo PBMC2sm > ../samples.txt
cat ../samples.txt

3a. Now, take a look at the raw data directory.

ls /share/workshop/intro_scrnaseq/$USER/scrnaseq_example/00-RawData

3b. To see a list of the contents of each directory.

ls *

3c. Lets get a better look at all the files in all of the directories.

ls -lah *

4. View the contents of the files using the ‘less’ command, when gzipped used ‘zless’ (which is just the ‘less’ command for gzipped files, q to exit):

Read 1

cd 00-RawData/
zless PBMC2sm_S18_L003_R1_001.fastq.gz
zless PBMC2sm_S18_L003_R2_001.fastq.gz

and Read 2

zless PBMC2sm_S18_L003_R2_001.fastq.gz

Make sure you can identify which lines correspond to a single read and which lines are the header, sequence, and quality values. Press ‘q’ to exit this screen. Then, let’s figure out the number of reads in this file. A simple way to do that is to count the number of lines and divide by 4 (because the record of each read uses 4 lines). In order to do this use cat to output the uncompressed file and pipe that to “wc” to count the number of lines:

zcat PBMC2sm_S18_L003_R1_001.fastq.gz | wc -l

Divide this number by 4 and you have the number of reads in this file. One more thing to try is to figure out the length of the reads without counting each nucleotide. First get the first 4 lines of the file (i.e. the first record):

zcat PBMC2sm_S18_L003_R1_001.fastq.gz  | head -2 | tail -1

Note the header lines (1st and 3rd line) and sequence and quality lines (2nd and 4th) in each 4-line fastq block. Then, copy and paste the DNA sequence line into the following command (replace [sequence] with the line):

echo -n [sequence] | wc -c

This will give you the length of the read.

Also can do the bash one liner:

echo -n $(zcat PBMC2sm_S18_L003_R1_001.fastq.gz  | head -2 | tail -1) | wc -c

See if you can figure out how this command works.

Quiz

Submit Quiz