1 RNA-seq原始数据

参考:https://mp.weixin.qq.com/s?__biz=MzA4NDAzODkzMA==&mid=2651272899&idx=1&sn=6779b2bd21f3b607a08227d31c7212c6&chksm=841ed2beb3695ba8bee23563c28caa005447b2298785719964732b16cafe3a15d7d4937b95c1&scene=21#wechat_redirect

原始数据格式为fastq,为本文文件.用来保存生物序列.每一个样本一个fastq文件,每个序列用四行信息记录.

2 数据整理

一般需要将不同数据放在不同文件夹下.

3 质量控制

使用fastQC和multiqc对测序质量进行评价.

fastqc -o <output dir> <seqfile1,seqfile2..>

fastqc -o qc_results *.fastq.gz

multiqc *fastqc.zip --pdf

4 使用trim-galore去除低质量的reads和adaptor

处理单个样本可以使用下面命令.

trim_galore -output_dir clean_data --paired --length 75 --quality 25 --stringency 5 /raw_data/*fastq.gz

批量多核处理:

首先设置路径到raw_data

ls|grep R1_001.fastq.gz>gz1
ls|grep R2_001.fastq.gz>gz2
paste gz1 gz2>config
cat config

然后创建一个trim.sh文件:

touch trim.sh

vim打开.

vim trim.sh

按esc,然后:i进入插入模式,写入下面的代码:

dir=/home/kelly/wesproject/clean_data/
cat config |while read id
do
      arr=${id}
      fq1=${arr[0]}
      fq2=${arr[1]}
      nohup trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done

其中dir是clean_data的绝对路径.

然后按esc,然后:wq保存退出.

然后运行该文件.

shell trim.sh

然后就会在clean_data中生成数据.

进入该文件夹下.可以查看文件的数据.

ls -lht | grep val | wc -l

这种写法|跟R中的%>%类似.

这里面,wc -l是计数的.

5 序列比对

我们使用hisat进行序列比对.

下载参考基因组:

路径设置到想要存放路径的地方.

wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg38.tar.gz

这里我们下载最近的hg38参考基因组.

hisat2使用.

hisat2 -p 6 -x <dir of index of genome> -1 seq_val_1.fq.gz -2 seq_val_2.fq.gz -S  tem.hisat2.sam

-p: 处理核心数.

-x: 参考基因组存放位置.最后需要写上geome.

-1: 两端测序的第一个文件.

-2: 两端测序的第二个文件.

-S: 生成的sam格式数据的名字.

具体例子:

hisat2 --dta -t -p 8 -x ./reference_genome/index/hg38/genome -1 ./clean_data/iPOP_MC_PBMC_RNAseq_1_S1_L001_R1_001_val_1.fq.gz -2 ./clean_data/iPOP_MC_PBMC_RNAseq_1_S1_L001_R2_001_val_2.fq.gz -S ./aligned/iPOP_MC_PBMC_RNAseq_1_S1_L001_R2_001_val_2.sam

路径需要在整个project root路径.

注意参考基因组的路径写法,其实我们的参考基因组就存放在./reference_genome/index/hg38/文件夹下.genome一定要加上去.

写一个循环进行批次处理:

for i in iPOP_MC_PBMC_RNAseq_1_S1_L001 iPOP_MC_PBMC_RNAseq_10_S10_L001 iPOP_MC_PBMC_RNAseq_11_S11_L001 iPOP_MC_PBMC_RNAseq_12_S12_L001 iPOP_MC_PBMC_RNAseq_13_S13_L001 iPOP_MC_PBMC_RNAseq_14_S14_L001 iPOP_MC_PBMC_RNAseq_15_S15_L001 iPOP_MC_PBMC_RNAseq_16_S16_L001 iPOP_MC_PBMC_RNAseq_17_S17_L001 iPOP_MC_PBMC_RNAseq_18_S18_L001 iPOP_MC_PBMC_RNAseq_19_S19_L001 iPOP_MC_PBMC_RNAseq_2_S2_L001 iPOP_MC_PBMC_RNAseq_20_S20_L001 iPOP_MC_PBMC_RNAseq_21_S21_L001 iPOP_MC_PBMC_RNAseq_22_S22_L001 iPOP_MC_PBMC_RNAseq_23_S23_L001 iPOP_MC_PBMC_RNAseq_24_S24_L001 iPOP_MC_PBMC_RNAseq_25_S25_L001 iPOP_MC_PBMC_RNAseq_26_S26_L001 iPOP_MC_PBMC_RNAseq_27_S27_L001 iPOP_MC_PBMC_RNAseq_28_S28_L001 iPOP_MC_PBMC_RNAseq_29_S29_L001 iPOP_MC_PBMC_RNAseq_3_S3_L001 iPOP_MC_PBMC_RNAseq_30_S30_L001 iPOP_MC_PBMC_RNAseq_31_S31_L001 iPOP_MC_PBMC_RNAseq_32_S32_L001 iPOP_MC_PBMC_RNAseq_33_S33_L001 iPOP_MC_PBMC_RNAseq_34_S34_L001 iPOP_MC_PBMC_RNAseq_35_S35_L001 iPOP_MC_PBMC_RNAseq_36_S36_L001 iPOP_MC_PBMC_RNAseq_37_S37_L001 iPOP_MC_PBMC_RNAseq_38_S38_L001 iPOP_MC_PBMC_RNAseq_39_S39_L001 iPOP_MC_PBMC_RNAseq_4_S4_L001 iPOP_MC_PBMC_RNAseq_40_S40_L001 iPOP_MC_PBMC_RNAseq_41_S41_L001 iPOP_MC_PBMC_RNAseq_42_S42_L001 iPOP_MC_PBMC_RNAseq_43_S43_L001 iPOP_MC_PBMC_RNAseq_44_S44_L001 iPOP_MC_PBMC_RNAseq_45_S45_L001 iPOP_MC_PBMC_RNAseq_46_S46_L001 iPOP_MC_PBMC_RNAseq_47_S47_L001 iPOP_MC_PBMC_RNAseq_5_S5_L001 iPOP_MC_PBMC_RNAseq_6_S6_L001 iPOP_MC_PBMC_RNAseq_7_S7_L001 iPOP_MC_PBMC_RNAseq_8_S8_L001 iPOP_MC_PBMC_RNAseq_9_S9_L001 Undetermined_S0_L001
do 
hisat2 --dta -t -p 15 -x ./reference_genome/index/hg38/genome \
-1 ./clean_data/"$i"_R1_001_val_1.fq.gz \
-2 ./clean_data/"$i"_R2_001_val_2.fq.gz \
-S ./aligned/"$i".sam; done

6 samtools将sam转换为bam

使用samtools将得到的sam格式数据转换为bam格式,并且进行sort.

单个转换:

samtools view -S iPOP_MC_PBMC_RNAseq_10_S10_L001.sam -b > iPOP_MC_PBMC_RNAseq_10_S10_L001.bam

批次转换:

for i in iPOP_MC_PBMC_RNAseq_11_S11_L001 iPOP_MC_PBMC_RNAseq_12_S12_L001 iPOP_MC_PBMC_RNAseq_13_S13_L001 iPOP_MC_PBMC_RNAseq_14_S14_L001 iPOP_MC_PBMC_RNAseq_15_S15_L001 iPOP_MC_PBMC_RNAseq_16_S16_L001 iPOP_MC_PBMC_RNAseq_17_S17_L001 iPOP_MC_PBMC_RNAseq_18_S18_L001 iPOP_MC_PBMC_RNAseq_19_S19_L001 iPOP_MC_PBMC_RNAseq_2_S2_L001 iPOP_MC_PBMC_RNAseq_20_S20_L001 iPOP_MC_PBMC_RNAseq_21_S21_L001 iPOP_MC_PBMC_RNAseq_22_S22_L001 iPOP_MC_PBMC_RNAseq_23_S23_L001 iPOP_MC_PBMC_RNAseq_24_S24_L001 iPOP_MC_PBMC_RNAseq_25_S25_L001 iPOP_MC_PBMC_RNAseq_26_S26_L001 iPOP_MC_PBMC_RNAseq_27_S27_L001 iPOP_MC_PBMC_RNAseq_28_S28_L001 iPOP_MC_PBMC_RNAseq_29_S29_L001 iPOP_MC_PBMC_RNAseq_3_S3_L001 iPOP_MC_PBMC_RNAseq_30_S30_L001 iPOP_MC_PBMC_RNAseq_31_S31_L001 iPOP_MC_PBMC_RNAseq_32_S32_L001 iPOP_MC_PBMC_RNAseq_33_S33_L001 iPOP_MC_PBMC_RNAseq_34_S34_L001 iPOP_MC_PBMC_RNAseq_35_S35_L001 iPOP_MC_PBMC_RNAseq_36_S36_L001 iPOP_MC_PBMC_RNAseq_37_S37_L001 iPOP_MC_PBMC_RNAseq_38_S38_L001 iPOP_MC_PBMC_RNAseq_39_S39_L001 iPOP_MC_PBMC_RNAseq_4_S4_L001 iPOP_MC_PBMC_RNAseq_40_S40_L001 iPOP_MC_PBMC_RNAseq_41_S41_L001 iPOP_MC_PBMC_RNAseq_42_S42_L001 iPOP_MC_PBMC_RNAseq_43_S43_L001 iPOP_MC_PBMC_RNAseq_44_S44_L001 iPOP_MC_PBMC_RNAseq_45_S45_L001 iPOP_MC_PBMC_RNAseq_46_S46_L001 iPOP_MC_PBMC_RNAseq_47_S47_L001 iPOP_MC_PBMC_RNAseq_5_S5_L001 iPOP_MC_PBMC_RNAseq_6_S6_L001 iPOP_MC_PBMC_RNAseq_7_S7_L001 iPOP_MC_PBMC_RNAseq_8_S8_L001 iPOP_MC_PBMC_RNAseq_9_S9_L001 Undetermined_S0_L001
do
samtools view -S "$i".sam -b > "$i".bam; done