转录组上游分析—使用iseq下载原始数据、小鼠基因组、单端测序数据处理
1 下载rawdata_使用iseq替代prefetch进行数据下载
进行数据集GSE105789上游分析的时候，总共才四个数据集，使用prefetch下载的时候，不知道网络抽了什么风，速度一直都很慢。下了10个小时才下了三分之一。!
经过大佬指点，尝试使用iseq替代prefetch进行数据下载
参考链接：https://mp.weixin.qq.com/s/AVqv07swFvjl6OCnLwwLPA
iseq可以直接下载GSE编号，不用再使用prefetch下载list中的SRA号，而且不需要再通过fastq-dump进行从SRA到fastq数据的转换。
首先我在当前conda环境中直接安装iseq是有冲突的，因此需要重新建立个conda环境：
正式使用
代码语言：bash复制conda create -n iseq
conda activate iseq
conda install iseq
#直接传输GSE编号即可，不用再使用prefetch下载list中的SRA号
iseq -i GSE105789 -g 
开始还是很快的，可以根据success.log中的内容确定哪些文件是已经成功下载的，
如果终止了重新下载需要把未下载完全的文件删除掉，如图中如果需要重新下载，需要删除掉SRR6203123,
同时可以看到下载下来的直接就是fastq.gz的文件，不必再进行sra to fastq的转换。
如图，重新下载会跳过已经下载好的。
后面的速度有所下降，可能是网络波动的问题，但是总体还是比prefetch快的。
这里还是采用prefetch下载(因为之前已经下载了90%)，然后通过fasterq-dump进行转换
单端数据
cat sra_data/SRR_download_List.txt | while read id
do
echo "fasterq-dump -e 32 --split-files -O fq_data/ --outfile ${id}.fastq sra_data/${id}/${id}.sra"
# 这一步使用fastq-dump或fasterq-dump都可以
echo "pigz -p 16 -f fq_data/${id}.fastq"
done > sra2fq.sh
less sra2fq.sh
nohup bash sra2fq.sh 1>sra2fq.log 2>&1 &
2 fastqc质控trim_galore过滤
代码语言：bash复制cd /home/data/t110558/project/GSE105789/data/rawdata/qc
nohup fastqc -t 20 -o ./ ../fq_data/SRR*.fastq.gz 1>qc.log 2>&1 &
#数据整合
multiqc *.zip -o ./ -n qc_fastqc 1>./multiqc.log 2>&1 &
#样本数量少就可以直接提交到后台运行
cd /home/data/t110558/project/GSE105789/data/rawdata/fq_data/
ls *gz|cut -d"." -f 1|sort -u |while read id;do
nohup trim_galore -j 20 -q 25 --phred33 --length 36 --stringency 3 --paired -o ../../cleandata ${id}*.gz 1>../../raw2clean.log 2>&1 &
done 
#修改代码
cd /home/data/t110558/project/GSE105789/data/rawdata/fq_data/
ls *gz|cut -d"." -f 1|sort -u |while read id;do
nohup trim_galore -j 8 -q 25 --phred33 --length 36 --stringency 3 -o ../../cleandata ${id}*.gz 1>../../raw2clean.log 2>&1 &
done 
3 hisat2比对
注意要使用小鼠的基因组
代码语言：bash复制cd /home/data/t110558/project/GSE105789/data/cleandata
indexPrefix=/home/data/t110558/database/GRCh38.111/Hisat2Index/GRCh38.dna
ls *gz|cut -d"_" -f 1|sort -u |while read id;do
nohup hisat2 -p 20 -x $indexPrefix -U ${id}*_trimmed.fq.gz -S ../../Mapping/hisat2/${id}.hisat.sam 1>../../Mapping/hisat2/hisat2.log 2>&1 &
done 
参考基因组选择小鼠
3.1 小鼠参考基因组下载
代码语言：bash复制# 下载基因组序列
nohup axel -n 20 https://ftp.ensembl.org/pub/release-112/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz >dna.log &
# 下载基因组注释文件
nohup axel -n 20 https://ftp.ensembl.org/pub/release-112/gtf/mus_musculus/Mus_musculus.GRCm39.112.gtf.gz >gtf.log &
nohup axel -n 20 https://ftp.ensembl.org/pub/release-112/gff3/mus_musculus/Mus_musculus.GRCm39.112.gff3.gz >gff.log&
#解压对应文件
nohup gunzip Mus_musculus.GRCm39.112.gtf.gz Mus_musculus.GRCm39.112.gff3.gz Mus_musculus.GRCm39.dna.primary_assembly.fa.gz >unzip.log &3.2 参考基因组建立索引
代码语言：bash复制cd /home/data/t110558/database/GRCm39.112
# Hisat2构建索引，构建索引时间比较长，建议提交后台运行，一般会运行20多分钟左右
vim Hisat2Index.sh
mkdir Hisat2Index
fa=Mus_musculus.GRCm39.dna.primary_assembly.fa
fa_baseName=GRCm39.dna
hisat2-build -p 50 -f ${fa} Hisat2Index/${fa_baseName}
# 运行
nohup bash Hisat2Index.sh >Hisat2Index.sh.log &
3.3 基因组比对
代码语言：bash复制cd /home/data/t110558/project/GSE105789/data/cleandata
indexPrefix=/home/data/t110558/database/GRCm39.112/Hisat2Index/GRCm39.dna
ls *gz|cut -d"_" -f 1|sort -u |while read id;do
nohup hisat2 -p 20 -x $indexPrefix -U ${id}*_trimmed.fq.gz -S ../../Mapping/hisat2/${id}.hisat.sam 1>../../Mapping/hisat2/hisat2.log 2>&1 &
done 
#sam转bam
ls *.sam|while read id ;do (nohup samtools sort -O bam -@ 20 -o $(basename ${id} ".sam").bam ${id} 1>sam2bam.log 2>&1 & );done
#为bam文件建立索引
ls *.bam |xargs -i samtools index {}
rm *.sam 
4 FeatureCounts定量
注意是单端数据
代码语言：bash复制cd /home/data/t110558/project/GSE105789/Expression/featureCounts
## 定义输入输出文件夹
gtf=/home/data/t110558/database/GRCm39.112/Mus_musculus.GRCm39.112.gtf
inputdir=/home/data/t110558/project/GSE105789/Mapping/hisat2
# featureCounts对bam文件进行计数
nohup featureCounts -T 20 -t exon -g gene_id -a $gtf -o all.id.txt $inputdir/*.bam 1>featureCounts.log 2>&1 &
# 对定量结果质控
multiqc all.id.txt.summary
# 得到表达矩阵
# 处理表头
#/home/data/t110558/project/GSE105789/Mapping/hisat2/SRR6203120.hisat.bam
less -S all.id.txt |grep -v '#' |cut -f 1,7- |sed 's#/home/data/t110558/project/GSE105789/Mapping/hisat2/##g' |sed 's#.hisat.bam##g' >raw_counts.txt
#查看表达矩阵
head -n 100 raw_counts.txt|column -t
