暴躁版
制约我们的是热情吗？是知识吗？
不，是网速啊，阿西吧！
1 获取SRA号——自己读文献，去pubmend搜
2 下载数据
cat SRR_Acc_List.txt | while read id; do (prefetch ${id} );done
别试了，龟速，下不下来，请愉快的打上一个&然后洗洗睡吧
jimmy贴心的给了一个kill命令的代码
ps -ef | grep prefetch | awk '{print $2}' | while read id; do kill ${id}; done
呵呵哒
3 质控
sra转fastq
for i in $wkd/*sra
do
echo $i
fastq-dump --split-3 --skip-technical --clip --gzip $i ## 批量转换
done
核心：fastq-dump
--split-3 确保一个双端测序的样本被拆分成两个fastq文件。
--skip-technical Dump only biological reads
-W|--clip Remove adapter sequences from reads
--gzip Compress output using gzip: deprecated, not
recommended （哈哈，为啥不推荐呢）
fastq转fastqc
ls *gz | xargs fastqc -t 4
multiqc ./
-t 线程数 （double哈哈）
还是好慢，rna-seq的正确打开方法是，找本书，一边seq一边看
multiqc一下，就得到了multiqc的报告，下下来看一下，打不开，loading report ，FU...K!马景涛咆哮版！
4.过滤
过滤质量差的reads和去除接头
ls ~/sra/*_1.fastq.gz >1
ls ~/sra/*_2.fastq.gz >2
paste 1 2 >config
cat 一下config,长这样
/trainee2/hz10/sra/SRR1039510_1.fastq.gz /trainee2/hz10/sra/SRR1039510_2.fastq.gz
/trainee2/hz10/sra/SRR1039511_1.fastq.gz /trainee2/hz10/sra/SRR1039511_2.fastq.gz
/trainee2/hz10/sra/SRR1039512_1.fastq.gz /trainee2/hz10/sra/SRR1039512_2.fastq.gz
trim_galore ,用于去除低质量和接头数据
改一下大神的脚本，换成自己的路径，长这样
conda activate rna2
bin_trim_galore=trim_galore
dir='/trainee2/hz10/sra/clean'
cat $1 |while read id
do
arr=(${id})
fq1=${arr[0]}
fq2=${arr[1]}
$bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o $dir $fq1 $fq2
done
conda deactivate rna2
然后bash一下，好慢，我忘了打&，然后掉线了，我恨！我需要把之前trim好的删掉再重新bash吗？还是它会很智能的自己跳过去呢？它好像跳过去了，nice~
太慢了，回家吃饭了。。。。。
好了，如下。
├── [ 123] 1
├── [ 123] 2
├── [ 246] config
├── [ 284] qc.sh
├── [2.9K] SRR1039510_1.fastq.gz_trimming_report.txt
├── [1.2G] SRR1039510_1_val_1.fq.gz
├── [3.1K] SRR1039510_2.fastq.gz_trimming_report.txt
├── [1.2G] SRR1039510_2_val_2.fq.gz
├── [2.9K] SRR1039511_1.fastq.gz_trimming_report.txt
├── [1.2G] SRR1039511_1_val_1.fq.gz
├── [3.1K] SRR1039511_2.fastq.gz_trimming_report.txt
├── [1.2G] SRR1039511_2_val_2.fq.gz
├── [2.9K] SRR1039512_1.fastq.gz_trimming_report.txt
├── [1.5G] SRR1039512_1_val_1.fq.gz
├── [3.1K] SRR1039512_2.fastq.gz_trimming_report.txt
└── [1.5G] SRR1039512_2_val_2.fq.gz
--phred33 Instructs Cutadapt to use ASCII+33 quality scores as Phred scores (Sanger/Illumina 1.9+ encoding) for quality trimming.
--stringency Overlap with adapter sequence required to trim a sequence. 剪切掉的overlap的序列
***--length <INT> **** Discard reads that became shorter than length INT because of either quality or adapter trimming. A value of '0' effectively disables this behaviour. Default: 20 bp.
--paired 双端测序 并且需要双端测序的两条序列都比length的设定值要长，可以再不影响fastq格式的情况下，过滤掉过短的序列。
4. 比对
使用hisat2
构建hisat2基因组索引,来自生信技能树论坛，传说很慢
# 构建目录
mkdir /mnt/d/reference/index
mkdir /mnt/d/reference/index/hisat
#下载
wget -P /mnt/d/reference/index/hisat [url=ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz]ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz[/url]
#解压，-C指定解压目录
tar -zxvf /mnt/d/reference/index/hisat/hg19.tar.gz -C /mnt/d/reference/index/hisat
# 移除下载的压缩包
rm /mnt/d/reference/index/hisat/*.tar.gz
#查看解压文件
ll /mnt/d/reference/index/hisat/hg19
Tips_1:不用单独给下载的hg19的index文件再分别构建目录，解压的时候会自动创建hg19目录。
Tips_2: 解压的文件中，包含genome.*.ht2的8个文件和一个shell脚本。
hisat2用法
hisat2 [options]* -x <hisat2-idx>{-1 <m1> -2 <m2> | -U <r> } [-S <hit>]
-x 指定基因组索引
-1 指定第一个fastq文件:双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。
-2 指定第二个fastq文件:双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。
-S 指定输出的SAM文件。
改了老师的批量运行代码
cd $wkd/clean
ls *gz|cut -d"_" -f 1 |sort -u |while read id;do
ls -lh ${id}_1_val_1.fq.gz ${id}_2_val_2.fq.gz
hisat2 -p 10 -x ~/sra/index/hisat/hg38/genome -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz -S ${id}.hisat.sam
done
运行了一下，运行完第一个不错，比对98%，第二个卡住了，报错，google了一下原因，sam文件太大，内存不够，服务器爆掉了。。。。我太难了。。。。
然后就一个一个做，运行完了转sam，或者直接转sam。
hisat2 -p 1 -x ~/sra/index/hisat2/hg38/genome -1 SRR1039511_1_val_1.fq.gz -2 SRR1039511_2_val_2.fq.gz | samtools sort -@ 1 -o ~/sra/clean/SRR1039511.sort.bam -
bam 转 sam
samtools sort -O bam -@ 5 -o SRR1039510.hisat.bam SRR1039510.hisat.sam
-O 输出文件格式
-@ 线程数
最后是要操作的文件名，上面用管道符的是最后的-
我想买个服务器，我太难了。。。。
终于结束了，为bam文件建立索引
ls *.bam |xargs -i samtools index {}
ls *.bam |while read id ;do ( samtools flagstat -@ 1 $id > $(basename ${id} ".bam").flagstat );done
samtools view SAM转换为二进制对应的BAM格式
samtools sort 比对排序
samtools index 构建索引。对排序文件进行索引之后，有利于快速提取基因组重叠区域的比对结果
samtools flagstats 给出BAM文件的比对结果
samtools常用命令详解
SRR1039512.sort.bam.bam.flagstat文件可以直接cat查看比对结果
5.获得表达矩阵
gtf="/teach/database/gtf/gencode.v29.annotation.gtf.gz"
featureCounts -T 5 -p -t exon -g gene_id -a $gtf -o all.id.txt ~/sra/clean/*.bam 1>counts.id.log 2>&1 &
featurecounts 据说最大的优点是快
-a 输入GTF/GFF基因组注释文件
-p 这个参数是针对paired-end数据。Check validity of paired-end distance when counting read pairs. Use -d and -D to set thresholds.
-F 指定-a注释文件的格式，默认是GTF
-g 从注释文件中提取Meta-features信息用于read count，默认是gene_id
-t 跟-g一样的意思，其是默认将exon作为一个feature
-o 输出文件
-T 多线程数
然后我们就得到表达矩阵 all.id.txt 啦，就可以用R进行下游分析啦啦啦啦（啦啦啦你妹），limma,DESeq2，想用什么用什么，哪里不会点哪里~
ps，我觉得jimmy大佬有一个地方写错了，align下面没有bam文件，应该是$wkd/clean/*bam~ 觉得自己超棒，有没有~
一天的时间，跑完了一个流程，但是，只跑了3个样本，服务器就over了，sad
明天再跟着洲更大神的视频跑一个流程，我是不是就算入门了呢？
（不交钱是我最后的倔强~）
附：jimmy原教程
网友评论